© M. Kerscher 2022

1. Introduction

Statistical approaches are often used to characterise the large-scale structure of the galaxy distribution, wherein it is assumed that the distribution of galaxies is a realisation of a point process (see e.g. Neyman & Scott 1958; Peebles 1980). Observations give us the positions of galaxies in space. From this set of points, we estimate the moments of the point process, specifically the two-point correlation function ξ(r) or the two-point density

$$\begin{array}{c}\hfill {\varrho }_2(\mathit{x},\mathit{y})={\varrho }^2(1+\xi (|\mathit{x}-\mathit{y}|)),\end{array}$$(1)

the probability of finding two galaxies at x and y, where 𝜚 is the number density. In a homogeneous and isotropic point process, ξ(r) only depends on the separation r = |x − y|. To determine ξ(r) from a galaxy catalogue within a finite domain W ⊂ ℝ³, we use estimators. In cosmology, estimators based on random point sets are most commonly used. These rely on the data–data DD, data–random DR, and random–random RR pair-counts (see below for the definition).

The two-point correlation function of the galaxy distribution is often used to constrain models of structure and galaxy formation and to estimate parameters of cosmological models. In current and upcoming galaxy samples, the positions of millions and up to billions of galaxies will be observed (Dawson et al. 2013, BOSS; Ross et al. 2020, eBOSS; Abbott et al. 2022, DES; Aghamousa et al. 2016, DESI; Alarcon et al. 2021, PAUS; Ivezić et al. 2019, LSST; Amendola et al. 2018, Euclid). Fast and reliable methods for calculating the two-point correlation function are needed. These large samples allow us to reduce the statistical error from cosmic variance, but for the error budget of the two-point correlation function we also have to control the systematic errors. One systematic contribution is the error from random sets used in the pair-counts DR and RR. As an example, consider the baryon acoustic oscillations (BAOs) which lead to a peak in the two-point correlation function of galaxies at a scale of about 100 Mpc h⁻¹ (Eisenstein et al. 2005; Bautista et al. 2021). This BAO peak has a height of approximately 0.01 above zero (compare Fig. 2 in Eisenstein et al. 2005). For a percent level accuracy, we need to calculate the two-point correlation function with an absolute accuracy of less than 10⁻⁴. Below, we show how to reduce the systematic error to this level of accuracy without sacrificing efficiency.

Several estimators for the two-point correlation function have been developed (Peebles & Hauser 1974; Hewett 1982; Davis & Peebles 1983; Rivolo 1986; Landy & Szalay 1993; Hamilton 1993). By comparing these estimators to a reference result from a cosmological simulation, Kerscher et al. (2000) found that the Landy & Szalay (1993) estimator is the preferred estimator with the smallest deviation from the reference and also that its bias is negligible compared to its variance.

We focus on methods for increasing the numerical accuracy of the pair-counts as used in these estimators. The random point set, shared in the pair-counts RR and DR is used to correct for boundary (finite-size) and inhomogeneous sampling effects. As we outline below, RR and DR are Monte Carlo volume integration schemes. As expected from standard Monte Carlo integration, the error of these pair-counts RR and DR scales at least as $1/\sqrt{N_{\text{r}}}$, where N_r is the number of random points used (see Sect. 3 for a more differentiated view). This slow convergence rate makes increasing the accuracy costly, and sometimes unfeasible. To improve the accuracy of the pair-counts without sacrificing efficiency, we follow two directions:

• The pair-counts can be expressed as averages of specific volume fractions. We use this to propose special adapted volume integration schemes which, in turn, can be calculated more efficiently than the standard approach.

• We replace the standard Monte Carlo scheme with a quasi-Monte Carlo integration, which leads to an improved scaling of the error that is almost proportional to 1/N_q, where N_q is the number of points from a low-discrepancy sequence.

We compare the standard and the new methods to exactly known reference values. This allows us to empirically validate the asymptotic scaling of the errors.

A variety of approaches have been suggested to improve the speed and accuracy of estimators for the two-point correlation function. Keihänen et al. (2019) show that at fixed computational cost, a split random catalogue improves the accuracy of estimators for the two-point correlation function. For galaxy catalogues, Demina et al. (2018) achieve a speedup by factorising the calculations in radial (redshift) and angular coordinates (see also Breton & de la Torre 2021). Perhaps closest to our work are the investigations by Dávila-Kurbán et al. (2021). These authors use glass-like point sets instead of the random catalogues, where we use low-discrepancy sequences. As reference values in our comparisons, we use the exact results from Baddeley et al. (1993) and Kerscher (1999) for a rectangular box as summarised in Appendix B (for periodic boxes see Appendix A.1). He (2021) also discusses some approximations for these exact results.

Other approaches focus on the computational problem of calculating the pair-counts. Tree-based methods can be significantly faster than a direct implementation of the pair-counts, specifically for small radii (Moore et al. 2001, see also Appendix C). The double loop in the pair-count calculations can be parallelised. Alonso (2012) showed how to obtain a speedup by a factor of 100 over the direct implementation by using multi-threading on multi-core CPUs or utilising many cores in GPUs. It is well known from matrix computations that the memory layout of the data can have dramatic consequences for the run time of algorithms (see e.g. Anderson et al. 1999). Also, for the pair-counts, a clever layout of the coordinates in the memory can lead to a significant speedup (Donoso 2019). A similar approach can be combined with multi-threading and vectorisation resulting in a blazingly fast code (Sinha & Garrison 2020). Our conceptual improvements can be combined with these computational speedups.

In Sect. 2 we give the definition of the pair-counts and discuss the geometry of the expected pair-counts. In Appendix A we give the details of the derivations and in Appendix B we summarize some results for simple sample geometries. Together with Sect. 2, this enables us to calculate exact reference values for the pair-counts. At the end of Sect. 2, we give a short introduction to the quasi-Monte Carlo method, as used in our improved estimators. In Sects. 3 and 4, we compare the standard and improved versions of the pair-counts with the exactly known reference values. We discuss the scaling of the error with the number of points. We put this together in Sect. 5 and show how an improved version of the Landy & Szalay (1993) estimator can be constructed. Again we discuss the scaling of the error with the number of points. In Sects. 5.1 and 5.2, we compare the run times of the new estimator with those of the standard Landy & Szalay (1993) estimator in some typical situations. In Sect. 6, we discuss how these improvements have to be adapted to estimate the two-point correlation function from an inhomogeneous sampled galaxy distribution. We summarise in Sect. 7 and give some recommendations. In Appendix C, we discuss details of the implementation, the run times, and give a link to the code.

2. Pair-counts, geometry, and quasi-Monte Carlo

The set of the N data points (e.g. galaxies) is ${\left\{{\mathit{x}}_i\right\}}_{i=1}^N$ with all points x_i ∈ W inside the observation window W ⊂ ℝ³, that is, inside the unmasked area. The number density is estimated with $\widehat{\rho }=\frac{N}{|W|}$, where |W| is the volume of W. We then define

$$\begin{array}{c}\hfill \mathsf{DD}(r)=\frac{1}{N^2}{\displaystyle \sum _{i=1}^N\sum _{j=1,j\ne i}^N\delta k_r^{\delta }(|{\mathit{x}}_i-{\mathit{x}}_j|),}\end{array}$$(2)

the normalised number of data–data pairs with a distance of r = |x_i − x_j| in the interval [r, r + δ]. We use a rectangular kernel

$$\begin{array}{c}\hfill k_r^{\delta }(s)=\frac{1}{\delta }{\mathbb{1}}_{[r,r+\delta ]}(s),\end{array}$$(3)

with the indicator function of the set A defined as

$$\begin{array}{c}\hfill {\mathbb{1}}_A(q)=\{\begin{array}{cc}1\hfill & \mathrm{if}q\in A,\hfill \\ 0\hfill & \mathrm{else}.\hfill \end{array}\end{array}$$(4)

Also, other kernels with $\int k_r^{\delta }(s)\text{d}s=1$ are possible (e.g. triangular, truncated Gaussian, or Epanechnikov). We consider N_r randomly distributed points ${\{{\mathit{y}}_j\}}_{j=1}^{N_{\text{r}}}$, all inside the sample geometry y_j ∈ W. The normalised number of data-random pairs with a distance in [r, r + δ] is denoted by

$$\begin{array}{c}\hfill \mathsf{DR}(r)=\frac{1}{N{N}_{\mathrm{r}}}{\displaystyle \sum _{i=1}^N\sum _{j=1}^{N_{\mathrm{r}}}\delta k_r^{\delta }(|{\mathit{x}}_i-{\mathit{y}}_j|).}\end{array}$$(5)

Similarly,

$$\begin{array}{c}\hfill \mathsf{RR}(r)=\frac{1}{{N_{\mathrm{r}}}^2}{\displaystyle \sum _{i=1}^{N_{\mathrm{r}}}\sum _{j=1,j\ne i}^{N_{\mathrm{r}}}\delta k_r^{\delta }(|{\mathit{y}}_i-{\mathit{y}}_j|),}\end{array}$$(6)

is the normalised number of random–random pairs. The Landy & Szalay (1993) estimator is defined as

$$\begin{array}{c}\hfill {\xi }_{\mathrm{LS}}(r)=\frac{\mathsf{DD}(r)-2\mathsf{DR}(r)+\mathsf{RR}(r)}{\mathsf{RR}(r)}·\end{array}$$(7)

Also, the estimators provided by Peebles & Hauser (1974), Hewett (1982), Davis & Peebles (1983), and Hamilton (1993) can be defined in terms of the pair-counts and our results apply accordingly.

2.1. Geometry of pair-counts

The expectation of the pair-counts DR and RR can be expressed in terms of geometric quantities depending on the sample window W and on the point set (for DR, Kerscher 1999). We first consider the set-covariance

$$\begin{array}{cc}\hfill {\gamma }_W(\mathit{x})& =|W\cap W_{\mathit{x}}|={\displaystyle {\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W(\mathit{y}){\mathbb{1}}_W(\mathit{y}+\mathit{x})\mathrm{d}\mathit{y},}\hfill \end{array}$$(8)

where W_x is the shifted window W, that is, the set of all points from W shifted by the vector x. |W ∩ W_x| is the volume of the set W ∩ W_x. The isotropised set-covariance $\overline{{\gamma }_W}(r)$ can be calculated from γ_W(x):

$$\begin{array}{c}\hfill \overline{{\gamma }_W}(r)=\frac{1}{4\pi }{\displaystyle {\int }_0^{\pi }{\int }_0^{2\pi }{\gamma }_W(\mathit{x}(r,\theta ,\varphi ))sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi .}\end{array}$$(9)

Here x(r, θ, ϕ) = (rcos(ϕ)sin(θ),rsin(ϕ)sin(θ),rcos(θ)). For a large number of random points N_r one obtains (see Appendix A):

$$\begin{array}{cc}\hfill \mathsf{RR}(r)\to \mathcal{RR}(r)& =\frac{4\pi }{{|W|}^2}{\displaystyle {\int }_r^{r+\delta }\overline{{\gamma }_W}(s)s^2\mathrm{d}s}\hfill \\ \hfill & \approx \frac{4\pi r^2\delta }{{|W|}^2}\overline{{\gamma }_W}(r)\mathrm{for}\delta \mathrm{small},\hfill \end{array}$$(10)

where ℛℛ(r) is the expectation value of the pair-count RR(r) illustrating its geometric nature.

Ripley (1988) used a local area weight in an estimator for his K-function (the normalised integrated two-point density) and Rivolo (1986) considered a similar weight in his estimator for the two-point correlation function. This weight is inversely proportional to the fraction of the surface area of a sphere B_r(y) with radius r centred on the point y inside W:

$$\begin{array}{c}\hfill \mathrm{area}(\partial B_r(\mathit{y})\cap W)={\displaystyle {\int }_0^{\pi }{\int }_0^{2\pi }{\mathbb{1}}_W(\mathit{y}+\mathit{x}(r,\theta ,\varphi ))sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi .}\end{array}$$(11)

For a large number of random points N_r, one obtains (see Appendix A):

$$\begin{array}{cc}\hfill \mathsf{DR}(r)\to \mathcal{DR}(r)& =\frac{1}{|W|N}{\displaystyle \sum _{i=1}^N{\int }_r^{r+\delta }\mathrm{area}(\partial B_s({\mathit{x}}_i)\cap W)\mathrm{d}s}\hfill \\ \hfill & \approx \frac{1}{|W|N}{\displaystyle \sum _{i=1}^N\mathrm{area}(\partial B_r({\mathit{x}}_i)\cap W)\delta ,\mathrm{for}\delta \mathrm{small}.}\hfill \end{array}$$(12)

As before, 𝒟ℛ(r) is the expectation value of the pair-count DR(r) illustrating its geometric nature. However, now both DR(r) and 𝒟ℛ(r) are depending on the points ${\left\{{\mathit{x}}_i\right\}}_{i=1}^N$ under consideration.

In Appendix B, we give expressions for $\overline{{\gamma }_W}(r)$ and area(∂B_r(y)∩W) if W is a rectangular box or a sphere. From these expressions, we calculate the reference values 𝒟ℛ and ℛℛ. This allows us to compare different integration schemes for the pair-counts DR and RR and we can investigate the scaling of the accuracy with the number of points used in these methods.

2.2. Quasi-Monte Carlo

In our improved methods for estimating the pair-counts, we use quasi-Monte Carlo integration. In a standard Monte Carlo integration scheme, one uses N_r random points {y₁, …, y_Nr} to estimate the integral ∫_[0, 1]^df(x)dx by $\frac{1}{N_{\text{r}}}\sum _{i=1}^{N_{\text{r}}}f({\mathit{y}}_i)$. The accuracy can be estimated using the Chebyshev inequality, which tells us that the probability of an error exceeding a given threshold decreases with $1/\sqrt{N_{\text{r}}}$. In other words, the standard error of a Monte Carlo integration scales as $1/\sqrt{N_{\text{r}}}$.

With quasi-Monte Carlo methods, we use $\frac{1}{N_{\text{q}}}\sum _{i=1}^{N_{\text{q}}}f({\mathit{q}}_i)$ to numerically integrate ∫_[0, 1]^df(x)dx (see e.g. Niederreiter 1992). This estimate of the integral almost looks identical to the Monte Carlo integration above. However, for a quasi-Monte Carlo integration, the points in Q = {q₁, …, q_Nq} are not random. It is essential for the application of quasi-Monte Carlo integration that for a given point set Q the error bound

$$\begin{array}{c}\hfill {\displaystyle |{\int }_{{[0,1]}^d}f(\mathit{x})\mathrm{d}\mathit{x}-\frac{1}{N_{\mathrm{q}}}\sum _{i=1}^{N_{\mathrm{q}}}f({\mathit{q}}_i)|\le V(f)D(Q)}\end{array}$$(13)

factorises into a measure of variation V(f) depending only on properties of f, and a measure of discrepancy D(Q) depending only on the properties of the point set Q. If we consider functions of bounded variation, Eq. (13) is referred to as the Koksama-Hlawka bound and the measure of discrepancy is the star discrepancy (see e.g. L’Ecuyer & Lemieux 2002). To control the error bound (Eq. (13)), we have to control D(Q) (we note that V(f) does not depend on the point set Q). Low-discrepancy sequences, such as the Halton sequence, have been constructed with that in mind. For such sequences, Halton (1960) showed that

$$\begin{array}{c}\hfill D(Q)\propto \frac{{(log{N}_{\mathrm{q}})}^d}{N_{\mathrm{q}}}·\end{array}$$(14)

For small dimensions, d, this compares favourably to a Monte Carlo integration where the standard error only scales as $1/\sqrt{N_{\text{r}}}$. Halton sequences can be used to estimate integrals over indicator functions, which in turn define volumes like the set covariance. As indicator functions have bounded total variation, Eq. (13) applies.

Upper bounds like Eq. (13) are worst-case bounds. Owen & Rudolf (2021) derive an analogue to a strong law of large numbers for randomised low-discrepancy sequences. This further justifies the procedure for scrambling the Halton sequence developed by Owen (2017), where the scaling from Eq. (14) in Eq. (13) still gives an upper bound, but on average smaller errors are expected. We use these randomised Halton sequences in our calculations (see Appendix C).

3. RR

First we investigate how the accuracy of the standard RR, as given in Eq. (6), scales with the number of random points used. The numerical implementation of the pair-counts is discussed in Appendix C. The expectation value ℛℛ of RR can be expressed in terms of the isotropised set-covariance; see Eq. (10). Using the Eqs. (B.1) and (B.3) for the isotropised set covariance, we calculate ℛℛ(r) as a reference value for rectangular boxes W.

Figure 1 shows the scaling of the relative error ε_RR = |ℛℛ(r)−RR(r)|/ℛℛ(r) with the number of random points used. The error ε_RR is the mean value calculated from 100 samples of the random point sets. Only for the largest N_r do we use 10 samples. This also applies to the errors calculated in the following sections. We show the results for a rectangular cuboid W = [0, 1]³ and for a thin rectangular box W = [0, 1]×[0, 0.1]². With the thin box, we investigate how the error is affected in a sample W where the boundary effects are more dominant. Landy & Szalay (1993, their Eq. (21)) show that the relative error of RR is of the form $\frac{A}{\sqrt{N_{\text{r}}}}+\frac{B}{N_{\text{r}}}$, with A and B depending on the window W and the bin width. As expected, we see in Fig. 1 that the error scales as $1/\sqrt{N_{\text{r}}}$ for large N_r. For small N_r and small radii, the contribution proportional to 1/N_r appears. This additional contribution can be seen in the voluminous sample but not in the thin box. Simply using a three-dimensional low-discrepancy sequence instead of the random points in Eq. (6) is not feasible, as we see below.

Fig. 1. Relative error εRR against the number of random points Nr used in the standard procedure for calculating RR(r) in rectangular windows W. The black dotted lines are proportional to 1/​Nr and $1/\sqrt{N_{\text{r}}}$.

From Eqs. (A.4) and (A.5) we get

$$\begin{array}{c}\hfill \mathcal{RR}(r)=\frac{1}{{|W|}^2}{\displaystyle {\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W(\mathit{y}){\mathbb{1}}_W(\mathit{z})\delta k_r^{\delta }(|\mathit{y}-\mathit{z}|)\mathrm{d}\mathit{y}\mathrm{d}\mathit{z}.}\end{array}$$(15)

This allows a more flexible approach. Consider two (random) point sets $P_1={\{{\mathit{y}}_i\}}_{i=1}^{N_{\text{r}}}$ and $P_2={\{{\mathit{z}}_j\}}_{j=1}^{N_{\text{r}}}$, with the points y_i ∈ W, z_j ∈ W all inside the sample geometry. We then define

$$\begin{array}{c}\hfill {\mathsf{RR}}_{2\mathrm{x}3\mathrm{d}}(r)=\frac{1}{{N_{\mathrm{r}}}^2}{\displaystyle \sum _{i=1}^{N_{\mathrm{r}}}\sum _{j=1}^{N_{\mathrm{r}}}\delta k_r^{\delta }(|{\mathit{y}}_i-{\mathit{z}}_j|),}\end{array}$$(16)

which is also an estimate of ℛℛ(r). If the points in P₁ and P₂ are drawn from a Poisson process, the estimate from Eq. (16) is almost the same as the estimate from Eq. (6), because points in a Poisson process are independent. See also Dávila-Kurbán et al. (2021), who observe that one cannot use P₁ ≡ P₂ for RR_2x3d when they construct an estimator for the two-point correlation function using ‘glass like’ point sets. The representation of ℛℛ in Eq. (15) as a six-dimensional integral suggests a six-dimensional (quasi-)Monte Carlo approach. Consequently, we use six-dimensional random points or a six-dimensional randomised Halton sequence ${\{{\mathit{q}}_i\}}_{i=1}^{N_{\text{q}}}$, which we split as q_i = (y_i, z_i) into two three-dimensional sequences. We scale the points in the sequences ${\{{\mathit{y}}_i\}}_{i=1}^{N_{\text{q}}}$ and ${\{{\mathit{z}}_i\}}_{i=1}^{N_{\text{q}}}$ such that each y_i ∈ W and z_j ∈ W are uniformly distributed inside the sample geometry.

Figure 2 compares the scaling of the relative error ε_RR, 2x3d = |ℛℛ(r)−RR_2x3d(r)|/ℛℛ(r) with the number of (quasi-)random points. For random points, we see the expected scaling $1/\sqrt{N}$ of the error. Using the quasi-Monte Carlo approach, we see that the error scales proportionally to 1/N, which is even faster than the theoretical expectation according to Eq. (14). Using a low-discrepancy sequence, we gain more than two orders of magnitude in accuracy compared to the random point sets.

Fig. 2. Relative error εRR, 2x3d calculated for rectangular windows W with a pure Monte Carlo integration (solid lines) and with a randomised Halton sequence (dashed lines). The black dotted lines are proportional to $1/\sqrt{N_{\text{r}}}$ (upper) and 1/Nq (lower).

4. DR

In Sect. 2.1 we show that the expectation of DR(r) is 𝒟ℛ(r) and that it can be calculated from the area fraction area(∂B_r(x_i)∩W). For a rectangular box, Baddeley et al. (1993) gave explicit expressions for area(∂B_r(x_i)∩W) (see Appendix B.1). We use these expressions to calculate 𝒟ℛ(r) according to Eqs. (A.9) and (A.10). In the numerical integration of Eq. (A.9), we make sure to achieve a relative error of at least 10⁻¹⁰ for our reference value 𝒟ℛ(r). Both DR(r) and 𝒟ℛ(r) depend on the point set under consideration. We need a realistic data set D = ${\left\{{\mathit{x}}_i\right\}}_{i=1}^N$ to calculate both the theoretical reference 𝒟ℛ(r) and the different estimates for DR(r). For this purpose, we use a sample of simulated galaxy clusters¹ from the Magneticum simulation (Hirschmann et al. 2014; Ragagnin et al. 2017). The side length of the simulation box is 325 Mpc h⁻¹ and we use the real space positions of the 10 429 simulated clusters. We rescale the coordinates by the side-length of the box such that all points are inside W = [0, 1]³. In the smaller window W = [0, 1]×[0, 0.1]² only 86 clusters are left.

Now we are set to determine the scaling of the relative error ε_DR = |𝒟ℛ(r)−DR(r)|/𝒟ℛ(r) with the number of (quasi-)random points used. Figure 3 compares the standard approach with ordinary random numbers in DR(r) to a pair-count determined using a low-discrepancy sequence instead of the random points. Using a randomised Halton sequence, we only observe a minor gain for small N, but for large N the error is reduced by an order of magnitude. In the small window, W = [0, 1]×[0, 0.1]², this is more pronounced. The scaling follows the expected behaviour from Eq. (14).

Fig. 3. Relative error εDR calculated for rectangular windows W with a standard Monte Carlo integration (solid lines) and with a quasi-Monte Carlo scheme using a randomised Halton sequence (dashed lines). The black dotted lines are proportional to $1/\sqrt{N}$ (upper) and (log N)3/N (lower).

We can do better if we consider the geometry of DR(r). From Eqs. (A.9) and (A.10) we know that the expectation value of DR(r) is $\mathcal{DR}(r)=\frac{1}{|W|N}\sum _{i=1}^N{\mathcal{V}}_r^{\delta }({\mathit{x}}_i)$, where

$$\begin{array}{cc}\hfill {\mathcal{V}}_r^{\delta }({\mathit{x}}_i)& ={\displaystyle {\int }_r^{r+\delta }\mathrm{area}(\partial B_s({\mathit{x}}_i)\cap W)\mathrm{d}s=|S_r^{\delta }({\mathit{x}}_i)\cap W|}\hfill \end{array}$$(17)

is the volume of the spherical shell with a radial range in [r, r + δ] around x_i inside the sample geometry W. As already suggested by Rivolo (1986), this directly leads to a Monte Carlo scheme. With N_sh points ${\{{\mathit{y}}_i\}}_{i=1}^{N_{\text{sh}}}$ (quasi-)randomly distributed in the shell $S_r^{\delta }({\mathit{x}}_i)=\{\mathit{y}\in {ℝ}^3|s<|\mathit{y}-{\mathit{x}}_i|\le s+\delta \}$ around x_i we define a (quasi-)Monte Carlo estimate of ${\mathcal{V}}_r^{\delta }({\mathit{x}}_i)$:

$$\begin{array}{cc}\hfill {\mathsf{V}}_s^{\delta }({\mathit{x}}_i)& =\frac{|S_r^{\delta }|}{N_{\mathrm{sh}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{sh}}}{\mathbb{1}}_W({\mathit{y}}_j),}\hfill \end{array}$$(18)

with the volume $|S_r^{\delta }|=\frac{4\pi }{3}({(r+\delta )}^3-r^3)$ of the shell. We then get ${\mathsf{V}}_s^{\delta }({\mathit{x}}_i)\to {\mathcal{V}}_s^{\delta }({\mathit{x}}_i)$ for a large number N_sh of (quasi-)random points. Consequently we compare

$$\begin{array}{c}\hfill {\mathsf{DR}}_{\mathrm{shell}}(r)=\frac{1}{|W|N}{\displaystyle \sum _{i=1}^N{\mathsf{V}}_r^{\delta }({\mathit{x}}_i)}\end{array}$$(19)

with 𝒟ℛ(r). This is not a pair-count, but we still have a double sum over N data points and now N_sh points in the shell.

Figure 4 shows ε_DR, shell = |𝒟ℛ(r)−DR_shell(r)|/𝒟ℛ(r). A comparison with DR in Fig. 3 shows that using DR_shell leads to a reduction of the error by almost two orders of magnitude even for ordinary random points. Again this can be improved by using low-discrepancy sequences; doing so allows us to additionally gain at least another order of magnitude.

Fig. 4. Relative error εDR, shell calculated for rectangular windows W with a pure Monte Carlo integration (solid lines) and with a randomised Halton sequence (dashed lines). The black dotted lines are proportional to $1/\sqrt{N_{\text{sh}}}$ (upper) and (log Nsh)3/Nsh (lower).

5. Estimating ξ

Now we join the improved pair-count estimates together and compare results from the standard Landy & Szalay (1993) estimator to our new estimator for the pair correlation function. We use the simulated galaxy clusters from the Magneticum simulation to illustrate the behaviour of the different estimators (Hirschmann et al. 2014; Ragagnin et al. 2017, compare also Sect. 4). The 10 429 clusters are in a box with side-length 325 Mpc h⁻¹ and we use their real space positions. The standard Landy & Szalay (1993) estimator was defined in Eq. (7):

$$\begin{array}{c}\hfill {\xi }_{\mathrm{LS}}(r)=\frac{\mathsf{DD}(r)-2\mathsf{DR}(r)+\mathsf{RR}(r)}{\mathsf{RR}(r)}·\end{array}$$

First we use the same N_r random points to calculate RR and DR (see Eqs. (5) and (6)). As an exact reference we have

$$\begin{array}{c}\hfill \mathrm{\Xi }(r)=\frac{\mathsf{DD}(r)-2\mathcal{DR}(r)+\mathcal{RR}(r)}{\mathcal{RR}(r)},\end{array}$$(20)

with the ℛℛ and 𝒟ℛ given in Eqs. (10), (12), using the results from Appendix B.1 for a rectangular box. In Sects. 3 and 4 we discuss alternative possibilities to calculate the pair-counts DR and RR. We focus on the following combination:

$$\begin{array}{c}\hfill {\stackrel{\sim }{\xi }}_{\mathrm{LS}}(r)=\frac{\mathsf{DD}(r)-2{\mathsf{DR}}_{\mathrm{shell}}(r)+{\mathsf{RR}}_{2\mathrm{x}3\mathrm{d}}(r)}{{\mathsf{RR}}_{2\mathrm{x}3\mathrm{d}}(r)},\end{array}$$(21)

which resembles the Landy & Szalay (1993) estimator, but now with improved pair-count estimates. We use a low-discrepancy sequence with N_shell 3D points to calculate DR_shell (see Eq. (19)), and another 6D low-discrepancy sequence with N_2x3d points to calculate RR_2x3d (see Eq. (16)).

As an illustrative example we compare these estimates in Fig. 5 using an insufficient number of points. As expected, we observe that a standard Landy & Szalay (1993) estimator with only 10⁴ random points shows deviations from the exact result, but with 10⁵ random points the LS-estimator starts to follow the exact result. Visually, one can see that the estimator using the low-discrepancy sequences achieves a higher accuracy already with N_2x3d = 10⁴ = N_sh points.

Fig. 5. Two-point correlation function calculated for a simulated galaxy cluster sample. The exact Ξ is compared to standard ξLS with Nr = 104 and Nr = 105 random points, and to the new ${\stackrel{\sim }{\xi }}_{\mathrm{LS}}$ using N2x3d = 104 = Nsh points from a randomised Halton sequence.

To make this quantitative and to disentangle the influence that these different methods for calculating DR and RR have, we consider the ‘half exact’ estimates for ξ(r):

$$\begin{array}{cc}\hfill {\mathrm{\Xi }}_{\mathcal{RR}}(r)& =\frac{\mathsf{DD}(r)-2\mathsf{DR}(r)+\mathcal{RR}(r)}{\mathcal{RR}(r)},\hfill \\ \hfill {\mathrm{\Xi }}_{\mathcal{DR}}(r)& =\frac{\mathsf{DD}(r)-2\mathcal{DR}(r)+\mathsf{RR}(r)}{\mathsf{RR}(r)},\hfill \end{array}$$(22)

where we use the exact reference values ℛℛ and 𝒟ℛ in turn. We also calculate Ξ_ℛℛ and Ξ_𝒟ℛ using DR_shell and RR_2x3d instead of DR and RR. For the comparison, we consider the absolute error Δ₁(r) = |Ξ(r)−Ξ_𝒟ℛ(r)| and Δ₂(r) = |Ξ(r)−Ξ_ℛℛ(r)|. In Fig. 6 we see that the scaling investigated in Sect. 3 directly transfers to the scaling of the error in the estimated ξ. For a fixed number of (quasi-)random points, RR_2x3d is more accurate. The gain in accuracy is considerable for intermediate and large radii. For small radii and for a smaller number of (quasi-)random points, the accuracy gain of RR_2x3d is reduced to a factor of two. With the standard DR, the error Δ₁ scales as $1/\sqrt{N_{\text{r}}}$ (see Fig. 7). Already using a low-discrepancy sequence instead of random points in DR leads to a reduction in the error and the scaling starts to follow 1/N_q. Compared to DR, the DR_shell gives a significantly smaller error, again showing a scaling proportional to 1/N_sh. This reduced error comes at a price: the run time of DR_shell is greater than the run time of DR for the same number of (quasi-)random points (see Appendix C).

Fig. 6. Absolute error Δ1 is shown for the standard RR (solid lines) and for RR2x3d (dashed lines) pair-counts. The black dotted lines are proportional to $1/\sqrt{N_{\text{r}}}$ (upper) and 1/N2x3d (lower).

5.1. Run time

Ultimately, we are not interested in the scaling of the error with the number of (quasi-)random points. We want to obtain reliable results fast. The speed of the calculation depends on the hardware platform, the algorithms, details of the implementation, and further effects (see Appendix C). Nevertheless, we prepared some examples and investigate what one can expect as a speedup when using the new estimator. In a given window, and for a given radius interval, we look at Figs. 6 and 7 to read off how many (quasi-)random points we approximately need to achieve a targeted absolute error. In most applications of the standard estimator ξ_LS, the same random point set is used with RR and DR. However, there is some room for optimisation. As we tune RR_2x3d and DR_shell in the new estimator separately, we also choose different sizes of the random samples in RR and DR to optimise the overall run time. With some provisional runs, we check that we achieve the targeted accuracy Δ within 10%; sometimes we had to adjust the numbers. Table 1 lists the final parameters we use. For case III, we use N_r = 10⁷ and achieve a Δ ≈ 5 × 10⁻⁴ for the standard estimator, and therefore the run times for case III may serve only as a lower bound.

Fig. 7. Absolute error Δ2 is shown for the standard DR (solid lines), the DR using a low-discrepancy sequence (dotted lines), and for DRshell also using a low-discrepancy sequence (dashed lines). The black dotted lines are proportional to $1/\sqrt{N_{\text{r}}}$ (upper) and 1/Nsh (lower).

The comparison of the run times was conducted on a small workstation (for details see Appendix C). We estimate the run times from 100 runs, for the longer running calculations we use 20 runs (in case III 5 runs). As the targeted accuracy is matched at the 10% level and this comparison did not run on an isolated machine, we can only give one digit for the run times. See Appendix C for further comments on the run time calculations.

The run times of the pair-counts are shown in Table 1 for different windows, radius ranges, and for different targeted accuracies. From Appendix C, we know that for a fixed number of (quasi-)random points, the standard pair-counts RR and DR are always faster than the corresponding RR_2x3d and DR_shell. However, we choose the number of (quasi-)random points by requiring a given accuracy. In that case, we need fewer points for the new pair-counts, which in turn lowers their run times. Now we compare RR and RR_2x3d as well as DR and DR_shell at fixed accuracy, and find that the run times of the new pair-counts are smaller in all the considered cases. We also see from Table 1 that the relative significance of the RR, RR_2x3d and DR, DR_shell parts depends on the considered case. For a targeted high accuracy (small Δ), the run times are dominated by the RR and RR_2x3d part, respectively. For small radii and for lower accuracy, the DR and DR_shell parts dominate the run time. The geometry of the window further differentiates this picture. In Table 2, we sum up the parts and show the run times of ξ_LS and the new ${\stackrel{\sim }{\xi }}_{\mathrm{LS}}$. In all the cases considered, we see a reduction in the run time using the new estimator. For high accuracy and for large and intermediate radii, the effect is considerable. Using the new estimator we expect a reduction in run time by a factor of between several hundred and almost 10⁴. For small radii and for less demanding accuracy goals, our method performs at least as well as the standard estimator and we may expect a reduction by up to an order of magnitude. The timing results for DD, RR, RR_2x3d and DR are independent of the number of radius intervals, because we can sort the counts into the radius bins. However, we need to repeat the DR_shell calculation for each radius bin.

5.2. A dense sample

Case VI from Table 1 shows that for small radii and lower accuracy goals, the run time is dominated by the DR or the DR_shell calculations. However, DR and DR_shell show very different run times depending on the number of data points (see Appendix C). Hence per se it is not clear which of the methods is preferable in dense samples for small radii and relaxed accuracy goals. We cannot check this with the cluster sample with N = 10429 data points only. Therefore, we use a sample of 1.36 million simulated galaxies² from the Magneticum simulation (Hirschmann et al. 2014; Ragagnin et al. 2017) in the box with side-length 325 Mpc h⁻¹. In Fig. 8, we see that the errors of RR and RR_2x3d both scale proportionally to 1/N_r (see the discussion in Sect. 3). For a fixed number of (quasi-)random points, the RR_2x3d is about twice as accurate as the RR for r ∈ [2, 3] Mpc h⁻¹ in this dense sample. However, the accuracies of DR and DR_shell differ drastically. Using the same number of points, accuracy is increased when using DR_shell instead of DR by at least a factor of 10³. We can understand this by looking at Eq. (19). In calculating DR_shell, we average the shell volume inside the window over all the points. With so many points in this dense sample, only a rough estimate of the shell volume for each individual point or galaxy is necessary to obtain a DR_shell of sufficient accuracy.

Fig. 8. Errors for estimates of the correlation function at r ∈ [2, 3] Mpc h−1 from 1.36 million galaxies and a subsample of 100 k galaxies in a box with side-length 325 Mpc h−1. Upper plot: absolute error Δ1 shown for the standard RR (solid lines) and for RR2x3d (dashed lines) pair-counts. Lower plot: absolute error Δ2 is shown for the standard DR (solid lines), the DR using a low-discrepancy sequence (dotted lines), and for DRshell also using a low-discrepancy sequence (dashed lines). The black dotted lines are proportional to $1/\sqrt{N_{\text{r}}}$ (upper) and 1/N2x3d (lower).

We proceed as in Sect. 5.1, choose a targeted accuracy Δ = 10⁻², and determine the required number of points N_r, N_2x3d, and N_sh in the pair-counts as summarised in Table 3. Extrapolating from Fig. 8, we could have guessed that the number N_sh of points, as used in DR_shell, is less than one. Instead of subsampling the galaxies, we simply use N_sh = 1, and therefore we are overestimating the required run time T_DR, shell. It is also interesting to observe that in case XI we only need N_r ≪ N random points for the pair-counts in the standard estimators to achieve the targeted accuracy.

From Table 3 we see that for the subsampled galaxy distribution (case X, 100 k points), the pair-counts RR and DR have a comparable contribution to the run time for the standard estimator, whereas for the new estimator, the RR_2x3d still provides the major contribution to the run time. For the full galaxy sample with 1.36 million galaxies, the dominant contribution is from DD followed by DR for the standard estimator. For the new estimator, both RR_2x3d and DR_shell contribute similarly to the run time. In both cases, the new estimator is still faster than the standard estimator, although only marginally in the full sample XI. These timing results should be considered with some reservation. As illustrated in Appendix C, a tree-based method can speedup the DD, RR, DR, and RR_2x3d for small radii significantly. For DR_shell, such a method for increasing the speed is currently unknown. This certainly affects the relative contribution of DR versus DR_shell and depending on the situation, DR using a tree-based method can be faster than DR_shell.

6. Incomplete sampling

The random points in RR and DR are not only used to correct for finite-size effects but also to correct for incomplete sampling. Here, we follow the ideas presented in Baddeley et al. (2000) and Shaw et al. (2021) and adapt them for the pair-counts. In this way, we can still use the new improved estimators and also correct for incomplete sampling.

Typically the galaxy distribution is observed incompletely. For example, unobserved regions around bright stars are masked. This leads to holes in the observational window W and can be dealt with using the methods already described. However, a partial sampling of the galaxy distribution in crowded fields leaves us with an inhomogeneous selection of the galaxy sample. Also, in magnitude-limited samples, we have a systematic selection of the galaxies, depending on the distance from us. We model this by p(x)∈(0, 1], the probability of including or observing a galaxy at position x. We assume that this p(x) is statistically independent of any other point in the galaxy distribution. For example, in a crowded field, only a fraction of the galaxies are randomly targeted for spectroscopic follow up. For a galaxy in this field, we have a p(x) equal to the fraction of targeted galaxies in the field. In a magnitude-limited sample, we observe all galaxies down to a limiting brightness l_lim. The luminosity of a galaxy L(d(x),l) at position x can be determined from the brightness l > l_lim and its luminosity distance d(x) from our galaxy (we neglect absorption for simplicity). Therefore, at a distance of d we only include galaxies with luminosity L > L(d, l_lim) in our catalogue. Now consider F(L), the distribution function of the absolute luminosities of all the galaxies (i.e. the normalised cumulative luminosity function). The fraction of galaxies included at a distance d(x) is then p(x) = 1 − F(L(d(x),l_lim)). We therefore need a good model for the luminosity function. The exclusion of a galaxy due to fibre collision cannot be modelled with the independent thinning (see below), and the following approach may only serve as an approximation. Further selection and sampling effects are discussed in Ross et al. (2012).

To investigate the effects of incomplete sampling, we assume an unobserved homogeneous and isotropic galaxy distribution G with number density 𝜚_G and two-point correlation function ξ_G(r). Our goal is to estimate this ξ_G(r). The inhomogeneous sampling, described by p(x), leads to an observed inhomogeneous galaxy catalogue and is modelled in the following way:

$$\begin{array}{c}\hfill {\mathsf{D}}^{\prime }=\{{\mathit{x}}_i\in G|{\mathit{x}}_i\in W\mathrm{and}{u}_i\le p({\mathit{x}}_i)\},\end{array}$$(23)

where u_i are independent random variables that are specifically independent of the points and uniformly distributed on [0, 1]. This closely follows the construction of an inhomogeneous Markov point process by independent thinning as discussed by Baddeley et al. (2000). The observed inhomogeneous point set ${\left\{{\mathit{x}}_i^{\prime }\right\}}_{i=1}^N$ is considered a realisation of D′. With the number density 𝜚_G of the homogeneous galaxy distribution, the inhomogeneous number density of this point process is 𝜚′(x) = p(x) 𝜚_G. The two-point density, the probability of observing a point at x and y, is

$$\begin{array}{c}\hfill {\varrho }_2^{\prime }(\mathit{x},\mathit{y})={\varrho }^{\prime }(\mathit{x}){\varrho }^{\prime }(\mathit{y})(1+{\xi }^{\prime }(|\mathit{x}-\mathit{y}|)).\end{array}$$(24)

This is well defined for a certain class of point process as defined in Baddeley et al. (2000). As the thinning is assumed to be independent from the points, we have

$$\begin{array}{c}\hfill {\xi }_G(|\mathit{x}-\mathit{y}|)={\xi }^{\prime }(|\mathit{x}-\mathit{y}|)=\frac{{\varrho }_2^{\prime }(\mathit{x},\mathit{y})}{{\varrho }_G^{2}p(\mathit{x})p(\mathit{y})}-1,\end{array}$$

for the two-point correlation function.

In full analogy to Sect. 2, we define the ‘inhomogeneous’ pair-counts DD′. The points ${\left\{{\mathit{x}}_i^{\prime }\right\}}_{i=1}^N$ are the incompletely sampled galaxies with positions ${\mathit{x}}_i^{\prime }$ inside the observation window W. Then

$$\begin{array}{c}\hfill {\mathsf{DD}}^{\prime }(r)=\frac{1}{N^2}{\displaystyle \sum _{i=1}^N\sum _{j=1,j\ne i}^N\delta k_r^{\delta }(|{\mathit{x}}_i^{\prime }-{\mathit{x}}_j^{\prime }|).}\end{array}$$(25)

To account for the effect of incomplete sampling, one often applies the thinning as described in Eq. (23) to the random points as well. In this way, our random point set ${\left\{{\mathit{y}}_i^{\prime }\right\}}_{i=1}^{N_{\text{r}}}$ is a realisation of an inhomogeneous Poisson process with number density $\frac{N_{\text{r}}}{N}{\varrho }^{\prime }(\mathit{x})$. The pair-counts involving these inhomogeneous random points are

$$\begin{array}{c}\hfill {\mathsf{RR}}^{\prime }(r)=\frac{1}{{N_{\mathrm{r}}}^2}{\displaystyle \sum _{i=1}^{N_{\mathrm{r}}}\sum _{j=1,j\ne i}^{N_{\mathrm{r}}}\delta k_r^{\delta }(|{\mathit{y}}_i^{\prime }-{\mathit{y}}_j^{\prime }|),}\end{array}$$(26)

and

$$\begin{array}{c}\hfill \mathsf{DR}\prime (r)=\frac{1}{N{N}_{\mathrm{r}}}{\displaystyle \sum _{i=1}^N\sum _{j=1}^{N_{\mathrm{r}}}\delta k_r^{\delta }(|{\mathit{x}}_i^{\prime }-{\mathit{y}}_j^{\prime }|).}\end{array}$$(27)

Appendix A.2 shows that the Peebles & Hauser (1974) estimator using DD′ and RR′ is ‘ratio unbiased’. For the Landy & Szalay (1993) estimator, we combine the pair-counts DD′, DR′, and RR′ with

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{LS}}^{\prime }(r)& =\frac{{\mathsf{DD}}^{\prime }(r)-2{\mathsf{DR}}^{\prime }(r)+{\mathsf{RR}}^{\prime }(r)}{{\mathsf{RR}}^{\prime }(r)},\hfill \end{array}$$(28)

in full analogy to Eq. (7). Now we show how we can calculate these inhomogeneous pair-counts using our new approach based on quasi-Monte Carlo methods. For this we make a detour and consider the expectation value of RR′ and then DR′.

For a large number of random points N_r from an inhomogeneous Poisson process, we obtain (see Appendix A.2):

$$\begin{array}{cc}\hfill {\mathsf{RR}}^{\prime }(r)& \to \mathcal{RR}\prime (r)=\frac{4\pi }{N^2}{\displaystyle {\int }_r^{r+\delta }\overline{{\mathrm{\Gamma }}_W}(s)s^2\mathrm{d}s,}\hfill \end{array}$$(29)

where $\overline{{\mathrm{\Gamma }}_W}(s)$ is given in Eq. (A.16), and ℛℛ′(r) is the expectation value of the pair-count RR′(r). There is no longer a direct correspondence with geometric objects like the set-covariance, however we can write ℛℛ′ as the following integral (compare to Eq. (15)):

$$\begin{array}{cc}\hfill \mathcal{RR}\prime (r)& =\frac{1}{N^2}{\displaystyle {\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W({\mathit{y}}_1){\mathbb{1}}_W({\mathit{y}}_2){\varrho }^{\prime }({\mathit{y}}_1){\varrho }^{\prime }({\mathit{y}}_2)}\hfill \\ \hfill & \times \delta k_r^{\delta }(|{\mathit{y}}_1-{\mathit{y}}_2|)\mathrm{d}{\mathit{y}}_1\mathrm{d}{\mathit{y}}_2.\hfill \end{array}$$(30)

This suggests a (quasi-)Monte Carlo approach similar to Eq. (16). We use homogeneously sampled point sets P₁, P₂ as described in Sect. 3. Then

$$\begin{array}{cc}\hfill {\mathsf{RR}}_{2\mathrm{x}3\mathrm{d}}^{\prime }(r)& =\frac{1}{N^2}{\displaystyle \sum _{i=1}^{N_{\mathrm{r}}}\sum _{j=1}^{N_{\mathrm{r}}}{\varrho }^{\prime }({\mathit{y}}_i){\varrho }^{\prime }({\mathit{z}}_j)\delta k_r^{\delta }(|{\mathit{y}}_i-{\mathit{z}}_j|),}\hfill \end{array}$$(31)

is an estimate of ℛℛ′(r) (we reiterate that the number density of the inhomogeneous Poisson process is $\frac{N_{\text{r}}}{N}{\varrho }^{\prime }(\mathit{x})$). As described in Sect. 3, the point sets P₁, P₂ can be generated from a six-dimensional random sequence or six-dimensional low-discrepancy sequence.

In Eq. (27) we use a thinned random point set to calculate DR′. Similarly, we could use a thinned low-discrepancy sequence in Eq. (27), but in Sect. 4 we see that this first approach is only mildly successful. We therefore follow the second approach by considering Eq. (A.18) and find

$$\begin{array}{cc}\hfill {\mathsf{DR}}^{\prime }(r)& \to \mathcal{DR}\prime (r)=\frac{1}{N^2}{\displaystyle \sum _{i=1}^N{\mathcal{U}}_r^{\delta }({\mathit{x}}_i),}\hfill \end{array}$$(32)

where ${\mathcal{U}}_r^{\delta }({\mathit{x}}_i)$ is given in Eq. (A.19). Similar to Eq. (18), we use N_sh points ${\{{\mathit{y}}_i\}}_{i=1}^{N_{\text{sh}}}$ (quasi-)randomly distributed in the shell $S_r^{\delta }({\mathit{x}}_i)=\{\mathit{y}\in {ℝ}^3|s<|\mathit{y}-{\mathit{x}}_i|\le s+\delta \}$ around x_i to estimate ${\mathcal{U}}_r^{\delta }({\mathit{x}}_i)$ by

$$\begin{array}{cc}\hfill {\mathsf{U}}_s^{\delta }({\mathit{x}}_i)& =\frac{|S_r^{\delta }|}{N_{\mathrm{sh}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{sh}}}{\mathbb{1}}_W({\mathit{y}}_j){\varrho }^{\prime }({\mathit{y}}_j).}\hfill \end{array}$$(33)

This results in a new estimate for 𝒟ℛ′(r):

$$\begin{array}{c}\hfill {\mathsf{DR}}_{\mathrm{shell}}^{\prime }(r)=\frac{1}{N^2}{\displaystyle \sum _{i=1}^N{\mathsf{U}}_r^{\delta }({\mathit{x}}_i).}\end{array}$$(34)

Similar to Eq. (21) we can use DD′, ${\mathsf{RR}}_{\text{2x3d}}^{\prime }$, and ${\mathsf{DR}}_{\text{shell}}^{\prime }$ to construct a Landy & Szalay (1993) – type estimator ${\stackrel{\sim }{\xi }}_{\mathrm{LS}}^{\prime }$ for inhomogeneously sampled galaxy catalogues.

Shaw et al. (2021) discuss how to estimate the ξ′(r) and a non-parametric model for 𝜚′(x) directly from the data. In a typical application to galaxy catalogues, we often have a good model for 𝜚′(x) (or equivalently for p(x)). Parameters of this model are fixed by the sampling strategy, and some are determined from the galaxy distribution (e.g. from an estimate of the luminosity function).

7. Summary and outlook

First, we focussed on the scaling of the error in estimates of the pair-counts RR and DR with the number of (quasi-)random points. For the standard approach, with ordinary random numbers, we confirm the expected slow shrinking of the error proportional to $1/\sqrt{N}$. A reformulation of the pair-counts makes a quasi-Monte Carlo integration possible. There we find that the error shrinks almost proportionally to 1/N_q, where N_q is the number of points from a low-discrepancy sequence. This scaling of the error not only holds for bulky samples but is even more pronounced in a thin sample with prevalent boundary effects. We are therefore confident that our improved methods are also applicable in more complicated sample geometries. We combine these improved pair-counts into a new Landy & Szalay (1993) – type estimator and compare with the standard one. The new estimator inherits the favourable scaling proportional to 1/N_q.

We can turn this observation around. For a fixed maximum error, we read off how many points are necessary to stay below this error with both the new and the standard estimator. We then compare the run times of the estimators. Depending on the accuracy goal, the radius range considered, and the density of the data, we observe a speedup by 50%, and up to a factor of almost 10⁴ for our improved estimator. More specifically, using RR_2x3d instead of the standard RR will increases the efficiency (accuracy vs. run time) of the estimator in any considered case. As tree-based methods are equally applicable to RR and RR_2x3d, we recommend using RR_2x3d for any use case. Also, in all the cases we considered, the DR_shell is faster than the standard DR at fixed accuracy. However, an unreserved recommendation is not possible. For the DR calculations, a tree based method could be employed to speed up the calculations (we did not use one in our calculations). Currently, no comparable method for a speedup is known for DR_shell. Therefore, in some situations for fixed accuracy, a tree-based method could lead to smaller run times in the DR calculations than in the DR_shell calculations. We expect this to be relevant for voluminous samples, small radii, and lower accuracy goals. A combination of DD, DR, RR_2x3d similar to Eq. (21) might therefore be useful. For large radii, we recommend DR_shell.

The random point sets are not only used for boundary corrections but to correct for incomplete sampling as well. Typically, the selection and sampling effects, as present in the galaxy distribution, are modelled onto the random point set used in the pair-counts. We discuss how to adapt this for the new improved estimator; we are using the probability of observing the galaxies as a weight in the calculation of the pair-counts. However, this is only a first step. Using weighting schemes, one could envision the construction of a minimum variance estimate for inhomogeneous sampled galaxies, similarly to Saunders et al. (1992), Feldman et al. (1994), and Colombi et al. (1998), but now for the improved estimators using low-discrepancy sequences.

Our use of the randomised Halton sequence is similar to the use of glass-like point sets by Dávila-Kurbán et al. (2021). Generating glass-like point sets can be computationally challenging. In contrast, a randomised Halton sequence is generated easily (see Appendix C). Halton sequences are probably the simplest choice, but other low-discrepancy sequences could be used as well (L’Ecuyer & Lemieux 2002). Interesting alternatives could be the so-called blue noise random sets, which are used in computer graphics for efficiently sampling from surfaces (see e.g. Heck et al. 2013). The computational improvements like tree-based methods or optimised memory layout are complimentary to our approach and can be used similarly for the pair-counts with low-discrepancy sequences. This should lead to a further speedup.

Acknowledgments

Many thanks to Adrian Baddeley for sharing code and the comments on the expressions for the area fraction. It is a pleasure to thank Klaus Dolag and Antonio Ragagnin for providing public access to the simulated galaxies and galaxy clusters from the Magneticum simulation. Also many thanks to the referee for his/her helpful comments and especially for the suggestion to include a more detailed discussion of the run times.

References

Appendix A: Expectation of pair-counts

Estimators for the two-point density and the correlation function using geometrical weights have been developed in spatial statistics (e.g. Stoyan & Stoyan 1994). These geometrical weights can be derived with an application of the Campbell-Mecke formula to the pair-counts. These ideas can be found in several places, such as for example Ripley (1976), Ohser (1983), Fiksel (1988), Stoyan et al. (1995) and quite recently for inhomogeneous point sets in Shaw et al. (2021). We provide analogous derivations to show the connection between the pair-counts DR and RR and geometric quantities such as the set-covariance and the area fraction (Kerscher 1999; Stoyan & Stoyan 2000).

The Campbell-Mecke formula connects the expectation value over realisations of a point process Φ to integrals over n-point densities (e.g. Stoyan et al. 1995; Daley & Vere-Jones 2003). For suitable functions f(x) and g(x, y), we have

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\displaystyle \sum _{\mathit{x}\in \mathrm{\Phi }}}f(\mathit{x})\right]& ={\displaystyle {\int }_{{\mathbb{R}}^3}f(\mathit{x})\varrho \mathrm{d}\mathit{x},}\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\displaystyle \sum _{\begin{array}{c}\mathit{x},\mathit{y}\in \mathrm{\Phi }\\ \mathit{x}\ne \mathit{y}\end{array}}}g(\mathit{x},\mathit{y})\right]& ={\displaystyle {\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}g(\mathit{x},\mathit{y}){\varrho }_2(\mathit{x},\mathit{y})\mathrm{d}\mathit{x}\mathrm{d}\mathit{y}.}\hfill \end{array}$$(A.2)

The Campbell-Mecke formula allows us to interchange the expectation of sums with an integration over the number density 𝜚 or the two-point density 𝜚₂(x, y). For a simple point process Φ, we define

$$\begin{array}{cc}\hfill S& ={\displaystyle \sum _{\begin{array}{c}\mathit{x},\mathit{y}\in \mathrm{\Phi }\\ \mathit{x}\ne \mathit{y}\end{array}}}{\mathbb{1}}_W(\mathit{x}){\mathbb{1}}_W(\mathit{y})\delta k_r^{\delta }(|\mathit{x}-\mathit{y}|).\hfill \end{array}$$(A.3)

Considering the galaxy distribution as a realisation of a point process Φ, we have $\mathsf{DD}=\frac{S}{N^2}$. We apply the Campbell-Mecke formula to 𝔼[S], where we assume homogeneity and isotropy 𝜚₂(x, y) = 𝜚₂(|x − y|).

$$\begin{array}{cc}\hfill \mathbb{E}[S]& ={\displaystyle {\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W(\mathit{x}){\mathbb{1}}_W(\mathit{y})\delta k_r^{\delta }(|\mathit{x}-\mathit{y}|){\varrho }_2(|\mathit{x}-\mathit{y}|)\mathrm{d}\mathit{x}\mathrm{d}\mathit{y}}\hfill \\ \hfill & ={\displaystyle {\int }_{{\mathbb{R}}^3}\underset{={\gamma }_W(\mathit{z})}{\underbrace{{\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W(\mathit{z}+\mathit{y}){\mathbb{1}}_W(\mathit{y})\mathrm{d}\mathit{y}}}\delta k_r^{\delta }(|\mathit{z}|){\varrho }_2(|\mathit{z}|)\mathrm{d}\mathit{z}}\hfill \\ \hfill & ={\displaystyle {\int }_0^{\infty }\underset{=4\pi \overline{{\gamma }_W}(s)}{\underbrace{{\int }_0^{\pi }{\int }_0^{2\pi }{\gamma }_W(\mathit{z}(s,\theta ,\varphi ))sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi }}\delta k_r^{\delta }(s){\varrho }_2(s)s^2\mathrm{d}s}\hfill \\ \hfill & =4\pi {\displaystyle {\int }_r^{r+\delta }\overline{{\gamma }_W}(s){\varrho }^2(1+\xi (s))s^2\mathrm{d}s,}\hfill \end{array}$$(A.4)

with 𝜚₂(s) = 𝜚²(1 + ξ(s)) and the number density 𝜚 of the point process.

The random point set used in RR is a realisation of a Poisson process with number density ${\varrho }_{\mathrm{r}}=\frac{N_{\text{r}}}{|W|}$ and by definition a vanishing two-point correlation function ξ(r) = 0. Applying the Campbell-Mecke formula to $\mathsf{RR}=\frac{S_{\mathrm{r}}}{{N_{\text{r}}}^2}$, we obtain the connection between the pair-count RR and the isotropised set-covariance (Kerscher 1999; Stoyan & Stoyan 2000):

$$\begin{array}{cc}\hfill \mathbb{E}[\mathsf{RR}(r)]=\mathcal{RR}(r)& =\frac{4\pi }{{N_{\mathrm{r}}}^2}{\displaystyle {\int }_r^{r+\delta }\overline{{\gamma }_W}(s){\varrho }_{\mathrm{r}}^2{s}^2\mathrm{d}s}\hfill \\ \hfill & \approx \frac{4\pi r^2\delta }{{|W|}^2}\overline{{\gamma }_W}(r)\mathrm{for}\delta \mathrm{small}.\hfill \end{array}$$(A.5)

Using the isotropised set covariance (see Appendix B.1), we can calculate ℛℛ(r) as a geometrical reference value for rectangular boxes.

As mentioned above, we consider the galaxy distribution as a realisation of a point process with number density $\varrho =\frac{N}{|W|}$ and we seek to estimate its two-point correlation function ξ(r). From

$$\begin{array}{cc}\hfill \frac{\mathbb{E}[\mathsf{DD}(r)]}{\mathbb{E}[\mathsf{RR}(r)]}& =\frac{{N_{\mathrm{r}}}^{2}4\pi {\int }_r^{r+\delta }\overline{{\gamma }_W}(s){\varrho }^2(1+\xi (s))s^2\mathrm{d}s}{N^{2}4\pi {\int }_r^{r+\delta }\overline{{\gamma }_W}(s){\varrho }_{\mathrm{r}}^2{s}^2\mathrm{d}s}\hfill \\ \hfill & =1+\frac{{\int }_r^{r+\delta }\overline{{\gamma }_W}(s)\xi (s)s^2\mathrm{d}s}{{\int }_r^{r+\delta }\overline{{\gamma }_W}(s)s^2\mathrm{d}s}\hfill \\ \hfill & \approx 1+\xi (r)\mathrm{for}\delta \mathrm{small},\hfill \end{array}$$(A.6)

we see that the estimator ${\xi }_{\mathrm{PH}}(r)=\frac{\mathsf{DD}(r)}{\mathsf{RR}(r)}-1$ of Peebles & Hauser (1974) is a ‘ratio unbiased’ estimator for the two-point correlation function ξ(r).

To derive an analogous relation between DR and the average surface area we consider

$$\begin{array}{cc}\hfill T& ={\displaystyle \sum _{i=1}^N\sum _{\mathit{y}\in \mathrm{\Phi }}{\mathbb{1}}_W(\mathit{y})\delta k_r^{\delta }(|{\mathit{x}}_i-\mathit{y}|),}\hfill \end{array}$$(A.7)

with a point process Φ and ${\left\{{\mathit{x}}_i\right\}}_{i=1}^N$ being a given set of points inside the sample geometry x_i ∈ W. With Φ, a Poisson process, and using the Campbell-Mecke formula for $\mathsf{DR}(r)=\frac{T}{N{N}_{\text{r}}}$, we obtain

$$\begin{array}{cc}\hfill \mathbb{E}[T]& ={\displaystyle \sum _{i=1}^N{\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W(\mathit{y})\delta k_r^{\delta }(|{\mathit{x}}_i-\mathit{y}|){\varrho }_{\mathrm{r}}\mathrm{d}\mathit{y}=\frac{N_{\mathrm{r}}}{|W|}\sum _{i=1}^N{\mathcal{V}}_r^{\delta }({\mathit{x}}_i),}\hfill \end{array}$$(A.8)

where ${\varrho }_{\mathrm{r}}=\frac{N_{\text{r}}}{|W|}$ is the number density of the Poisson process and

$$\begin{array}{cc}\hfill {\mathcal{V}}_r^{\delta }({\mathit{x}}_i)& =|S_r^{\delta }({\mathit{x}}_i)\cap W|={\displaystyle {\int }_r^{r+\delta }\mathrm{area}(\partial B_s({\mathit{x}}_i)\cap W)\mathrm{d}s}\hfill \end{array}$$(A.9)

is the volume of the spherical shell $S_r^{\delta }({\mathit{x}}_i)=\{\mathit{y}\in {ℝ}^3|s<|\mathit{y}-{\mathit{x}}_i|\le s+\delta \}$ with a radial range in [r, r + δ] around x_i inside the sample geometry W. For simple sample geometries, the area(∂B_s(x_i)∩W) can be calculated explicitly (see Appendix B.1). The integral in Eq. (A.9) can be evaluated using standard numerical methods to obtain the expectation value of value DR(r):

$$\begin{array}{cc}\hfill \mathcal{DR}(r)& =\frac{1}{|W|N}{\displaystyle \sum _{i=1}^N{\mathcal{V}}_r^{\delta }({\mathit{x}}_i)}\hfill \\ \hfill & \approx \frac{1}{|W|}\frac{1}{N}{\displaystyle \sum _{i=1}^N\mathrm{area}(\partial B_r({\mathit{x}}_i)\cap W)\delta \mathrm{for}\delta \mathrm{small}.}\hfill \end{array}$$(A.10)

Clearly, the expectation value of the estimators from Davis & Peebles (1983), Hewett (1982), and Hamilton (1993) can be expressed in terms of the istotropised set-covariance and the average area fraction (Kerscher 1999).

A.1. Periodic boundaries

Let us assume that our window is a rectangular box W = [0, L₁]×[0, L₂]×[0, L₃] with periodic boundaries, i.e. W has the topology of a three-torus. The majority of cosmological simulations enforce these boundary conditions. In such a situation, no boundary corrections are needed in the calculation of the two-point correlation function, but we have to respect the periodicity in each coordinate direction (Stoyan et al. 1995). The distance between two points x, y ∈ W is $d(\mathit{x},\mathit{y})=\sqrt{d_1^2+d_2^2+d_3^2}$, with d_i = min{|x_i − y_i|,L_i − |x_i − y_i|}. This is applicable if r is smaller than any L_i/2. For the shifted periodic box, we have W_x = {y | y − x ∈ W}=W and consequently the (isotropised) set-covariance is constant ${\gamma }_W(\mathit{x})=\overline{{\gamma }_W}(|\mathit{x}|)=|W|$. Similar to the derivation in Eq. (A.4), the expectation value of DD(r) can be calculated using the Campbell-Mecke formula:

$$\begin{array}{cc}\hfill \mathbb{E}[\mathsf{DD}(r)]& =\frac{4\pi }{N^2}{\displaystyle {\int }_r^{r+\delta }|W|{\varrho }^2(1+\xi (s))s^2\mathrm{d}s}\hfill \\ \hfill & =\frac{|S_r^{\delta }|}{|W|}+\frac{4\pi }{|W|}{\displaystyle {\int }_r^{r+\delta }\xi (s)s^2\mathrm{d}s}\hfill \\ \hfill & \approx \frac{4\pi r^2\delta }{|W|}(1+\xi (s))\mathrm{for}\delta \mathrm{small},\hfill \end{array}$$(A.11)

where the volume of the shell $|S_r^{\delta }|=\frac{4\pi }{3}({(r+\delta )}^3-r^3)$. Therefore, for small δ, $\frac{|W|}{|S_r^{\delta }|}\mathsf{DD}(r)-1$ is an unbiased estimate of the two-point correlation function ξ(r) in a periodic box; neither a random points set, nor a boundary correction with geometric factors is needed.

A.2. Expectation of pair-counts for inhomogeneous point sets

As described in Sect. 6 we are naturally confronted with an inhomogeneous sampled galaxy distribution. Here we calculate the expectation values of DD′,DR′ and RR′ for such an inhomogeneous situation.

One also has Campbell-Mecke formulas for an inhomogeneous point process Φ′ (Baddeley et al. 2000; Shaw et al. 2021):

$$\begin{array}{cc}\hfill \mathbb{E}{\displaystyle \left[\sum _{\mathit{x}\in \mathrm{\Phi }\prime }f(\mathit{x})\right]}& ={\displaystyle {\int }_{{\mathbb{R}}^3}f(\mathit{x}){\varrho }^{\prime }(\mathit{x})\mathrm{d}\mathit{x},}\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill \mathbb{E}{\displaystyle \left[\sum _{\begin{array}{c}\mathit{x},\mathit{y}\in \mathrm{\Phi }\prime \\ \mathit{x}\ne \mathit{y}\end{array}}g(\mathit{x},\mathit{y})\right]}& ={\displaystyle {\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}g(\mathit{x},\mathit{y}){\varrho }_2^{\prime }(\mathit{x},\mathit{y})\mathrm{d}\mathit{x}\mathrm{d}\mathit{y}}\hfill \\ \hfill & ={\displaystyle {\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}g(\mathit{x},\mathit{y}){\varrho }^{\prime }(\mathit{x}){\varrho }^{\prime }(\mathit{y})(1+{\xi }^{\prime }(|\mathit{x}-\mathit{y}|))\mathrm{d}\mathit{x}\mathrm{d}\mathit{y}.}\hfill \end{array}$$(A.13)

Again we consider

$$\begin{array}{cc}\hfill S^{\prime }& ={\displaystyle \sum _{\begin{array}{c}\mathit{x},\mathit{y}\in \mathrm{\Phi }\prime \\ \mathit{x}\ne \mathit{y}\end{array}}}{\mathbb{1}}_W(\mathit{x}){\mathbb{1}}_W(\mathit{y})\delta k_r^{\delta }(|\mathit{x}-\mathit{y}|),\hfill \end{array}$$(A.14)

and calculate the expectation (see Baddeley et al. 2000; Shaw et al. 2021 for similar derivations)

$$\begin{array}{cc}\hfill \mathbb{E}[S^{\prime }]& ={\displaystyle {\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W(\mathit{x}){\mathbb{1}}_W(\mathit{y})\delta k_r^{\delta }(|\mathit{x}-\mathit{y}|){\varrho }_2\prime (|\mathit{x}-\mathit{y}|)\mathrm{d}\mathit{x}\mathrm{d}\mathit{y}}\hfill \\ \hfill & ={\displaystyle {\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W(\mathit{z}+\mathit{y}){\mathbb{1}}_W(\mathit{y}){\varrho }^{\prime }(\mathit{z}+\mathit{y}){\varrho }^{\prime }(\mathit{y})\mathrm{d}\mathit{y}\delta k_r^{\delta }(|\mathit{z}|)(1+{\xi }^{\prime }(|\mathit{z}|)\mathrm{d}\mathit{z}}\hfill \\ \hfill & =4\pi {\displaystyle {\int }_r^{r+\delta }{\overline{\mathrm{\Gamma }}}_W(s)(1+{\xi }^{\prime }(s))s^2\mathrm{d}s,}\hfill \end{array}$$(A.15)

where ${\overline{\mathrm{\Gamma }}}_W(s)$ is the density weighted isotropised set-covariance

$$\begin{array}{cc}\hfill {\overline{\mathrm{\Gamma }}}_W(s)& =\frac{1}{4\pi }{\displaystyle {\int }_0^{\pi }{\int }_0^{2\pi }{\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W(\mathit{z}(s,\theta ,\varphi )+\mathit{y}){\mathbb{1}}_W(\mathit{y})}\hfill \\ \hfill & {\varrho }^{\prime }(\mathit{z}(s,\theta ,\varphi )+\mathit{y}){\varrho }^{\prime }(\mathit{y})\mathrm{d}\mathit{y}sin(\theta )\mathrm{d}\theta \mathrm{d}\varphi .\hfill \end{array}$$(A.16)

For a homogeneous point distribution with 𝜚′(x) = 𝜚, we consistently get ${\overline{\mathrm{\Gamma }}}_W(s)={\varrho }^2\overline{{\gamma }_W}(s)$.

The random point set used in RR′ is a realisation of an inhomogeneous Poisson process with number density $\frac{N_{\text{r}}}{N}{\varrho }^{\prime }(\mathit{x})$, and a vanishing two-point correlation function ξ′(r) = 0. N is the number of galaxies and N_r the number of random points in W. From Eq. (A.15) we get

$$\begin{array}{cc}\hfill \mathbb{E}[{\mathsf{RR}}^{\prime }(r)]=\mathcal{RR}\prime (r)& =\frac{4\pi }{N^2}{\displaystyle {\int }_r^{r+\delta }\overline{{\mathrm{\Gamma }}_W}(s)s^2\mathrm{d}s}\hfill \\ \hfill & \approx \frac{4\pi r^2\delta }{N^2}\overline{{\mathrm{\Gamma }}_W}(r)\mathrm{for}\delta \mathrm{small}.\hfill \end{array}$$(A.17)

Using Eq. (A.15) and Eq. (A.17) we can determine

$$\begin{array}{cc}\hfill \frac{\mathbb{E}[\mathsf{DD}\prime (r)]}{\mathbb{E}[{\mathsf{RR}}^{\prime }(r)]}& =\frac{{\int }_r^{r+\delta }\overline{{\mathrm{\Gamma }}_W}(s)(1+{\xi }^{\prime }(s))s^2\mathrm{d}s}{{\int }_r^{r+\delta }\overline{{\mathrm{\Gamma }}_W}(s)s^2\mathrm{d}s}\hfill \\ \hfill & \approx 1+{\xi }^{\prime }(r)\mathrm{for}\delta \mathrm{small},\hfill \end{array}$$

Hence we arrive at the well-known result that the estimator ${\xi }_{\mathrm{PH}}^{\prime }(r)=\frac{\mathsf{DD}\prime (r)}{{\mathsf{RR}}^{\prime }(r)}-1$ of Peebles & Hauser (1974) is ‘ratio unbiased’. ${\xi }_{\text{PH}}^{\prime }(r)$ is an estimate of the two-point correlation function ξ′(r) = ξ_G(r) if we apply the same selection effects to the random point set as we find in the incompletely sampled galaxy distribution.

For a given point set ${\left\{{\mathit{x}}_i\right\}}_{i=1}^N$ inside the sample geometry x_i ∈ W, we calculate the expectation of DR′. Using the Campbell-Mecke formula and with Φ′ an inhomogeneous Poisson process with number density $\frac{N_{\text{r}}}{N}{\varrho }^{\prime }(\mathit{x})$, we obtain

$$\begin{array}{cc}\hfill \mathbb{E}[{\mathsf{DR}}^{\prime }]& =\mathcal{DR}\prime (r)=\frac{1}{N{N}_{\mathrm{r}}}{\displaystyle \sum _{i=1}^N{\int }_{{\mathbb{R}}^3}{\mathbb{1}}_W(\mathit{y})\delta k_r^{\delta }(|{\mathit{x}}_i-\mathit{y}|)\frac{N_{\mathrm{r}}}{N}{\varrho }^{\prime }(\mathit{y})\mathrm{d}\mathit{y}}\hfill \\ \hfill & =\frac{1}{N^2}{\displaystyle \sum _{i=1}^N{\mathcal{U}}_r^{\delta }({\mathit{x}}_i),}\hfill \end{array}$$(A.18)