Corrfunc --- A Suite of Blazing Fast Correlation Functions on the CPU

1 INTRODUCTION

• The paper focuses on the 2PCF and outlines implementation, CPU hardware, optimizations, performance, comparisons, shortcomings, and future directions.Source: §1 · pp. 1–2
• The 2PCF probes both cosmology and galaxy formation, and larger galaxy samples can improve the precision of correlation functions.Upcoming surveys target tens of millions to billions of galaxies.Source: §1 · pp. 1–2
• Brute-force 2PCF computation scales as O(N^2), requiring 10^12 pair separations for 10^6 galaxies.For tens of millions of galaxies, the brute-force calculation can take days to weeks.Source: §1 · p. 2
• Repeated brute-force calculations for MCMC modelling could take at least a year for 10^5 iterations, making such analyses unfeasible at the stated timescale.The estimate assumes five minutes per 2PCF calculation.Source: §1 · p. 2
• Corrfunc aims to fill the gap between HPC codes that often omit astronomy use cases and astronomy codes described as slower and problem-specific.The paper presents it as a high-performance, well-tested, documented, flexible, open-source code for many correlation functions.Source: §1 · p. 2

2 BACKGROUND

• A 2PCF measures excess pair probability at spatial or angular separation, with spatial ξ(r) and angular ω(θ) forms.The definitions use mean galaxy density and differential volume or solid angle.Source: §2 · p. 2
• The projected correlation function wp(rp) integrates ξ(rp, π) along the line of sight to account for peculiar-velocity effects, using a finite πmax in practice.The passage gives a typical πmax range of 40–80 Mpc and notes the boundary trade-off.Source: §2 · p. 3

2.1 How to compute a correlation function

• Pair counting computes separations and bins them; Corrfunc’s design separates this work into distance, bin-index, and optional pair-weight functions.The natural estimator is described as 1 + ξ(r) = DD(r)/RR(r), while the Landy-Szalay estimator also uses a DR term [Landy & Szalay 1993].Source: §2.1 · p. 3

2.2 Package Summary

• Corrfunc emphasizes correctness, shared-memory OpenMP performance, portability, flexible access through APIs or command-line tools, and convenient installation.The paper describes a default installation requiring no user input while allowing advanced customization.Source: §2.2 · pp. 3–4
• Corrfunc provides routines for Cartesian simulated data and spherical-coordinate observed data, with different line-of-sight conventions for theory and mocks routines.The passage describes four simulated-galaxy routines and two observed-galaxy routines.Source: §2.2 · p. 4

3.1 Evolution of CPU design towards multi-core processors

• CPU clock-speed growth stalled as power dissipation became a constraint, shifting hardware toward multiple lower-frequency cores for increased throughput.The passage contrasts two cores with running one core at twice the frequency, which would require four times the power.Source: §3.1 · p. 4

3.2 The CPU-Memory Performance Gap: The Emergence of the Cache Hierarchy

• Cache hierarchies place progressively larger, slower memory closer to the CPU, and cache misses can substantially increase execution time.The example reports that a 5% lower cache hit rate can produce about a threefold runtime increase when misses require RAM access.Source: §3.2 · p. 5

3.3 Singe Instruction Multiple Data – SIMD (Vectorisation)

• SIMD executes independent array operations on blocks of values, and AVX can operate on eight floats simultaneously, potentially yielding an 8× speedup over un-vectorized code.The computation must be reformulated into independent chunks to realize such gains.Source: §3.3 · p. 5

3.3.1 How to Create Vectorised Software

• Corrfunc uses C with explicit SIMD intrinsics because histogram updates create dependencies that prevent automatic vectorization.The paper says SIMD libraries would have required additional C++ knowledge and delayed development.Source: §3.3.1 · p. 6

4 METHODS OVERVIEW

• Corrfunc combines reduced computation with predictable contiguous memory access, vectorized operations, and multicore parallelism.Its strategy prunes potential search volume before seeking further reductions at the particle level.Source: §4 · p. 6

4.1 Reducing the Total Number of Distance Computations: Partitioning the particles on a grid

• Fixed-width grid partitioning reduces the regions searched for pairs within Rmax and has lower partitioning overhead than the tested tree strategy.The paper reports that axis-aligned subdivisions outperformed the tree algorithms it initially tested.Source: §4.1 · pp. 6–7

4.1.1 Partitioning the particles on a grid: Spatial

• For spatial correlations, cells at least Rmax wide let the algorithm prune cells more than one cell away along any dimension.The method builds a bounding cuboid around both datasets; with equal grid dimensions, the described 2-D search covers nine cells and the 3-D search 27.Source: §4.1.1 · pp. 7–8, Figure 1 · p. 7

4.1.2 Partitioning the particles on a grid: Angular

• For angular correlations, Corrfunc bins declination and can further prune candidates with right-ascension bins whose allowed widths depend on declination.Near the poles it uses the entire RA range in a DEC bin, creating only one RA cell there.Source: §4.1.2 · p. 8, Figure 2 · p. 8

4.2 Reducing the Total Number of Distance Computations: Sub-dividing the Cells

• 84% of total separations computed in practice must be discarded with default-width cells; for 2-D Cartesian positions, only about 24% satisfy the distance cuts.The passage attributes the 84% figure to integer cell counts making cell-width exceed Rmax.Source: §4.2 · p. 9
• Subdividing each Rmax cell into k sub-cells reduces distance computations to (k+1)/2k of the original, approaching one-half, but increases neighbor-cell checks.In 3-D, the neighbor-cell count rises from 27 at (1,1,1) to 125 at (2,2,2), so k must balance computation savings against loop overhead.Source: §4.2 · p. 9

4.3 Reducing the Total Number of Distance Computations: Associating Pairs of Cells

• Spatial grids associate cells that can contain pairs within the maximum separation; periodic boundaries wrap target indices and store per-axis offsets for later separation calculations.With spatial cell widths at least the maximum separation, candidate pairs span at most one cell along each Cartesian axis.Source: §4.3 · pp. 9–10
• Angular grids with staggered right-ascension bands add clockwise and counter-clockwise cross-band cell pairs, then retain only unique target cells.Projected right-ascension boundaries identify the source cell's extent on the target declination band.Source: §4.3 · p. 10

4.4 Reducing the Total Number of Distance Computations: Symmetry considerations for Autocorrelations

• Auto-correlations count unique cell and particle pairs, double the resulting bin counts, and correct self-pairs when Rmin is zero.This makes auto-correlation results match cross-correlation of a dataset with itself, with an explicit correction to avoid counting zero-separation self-pairs twice.Source: §4.4 · p. 10

4.5 Speeding up the Calculation of Separations in Cell-Pairs: Improving Cache Locality

• Corrfunc duplicates particle data into contiguous per-cell arrays, increasing memory use to improve cache access during cell-pair processing.The passage estimates typical particle data at about 20–50 MB and says contiguous locations provide a significant performance boost on modern CPUs.Source: §4.5 · p. 10

4.6 Speeding up the Calculation of Separations in Cell-Pairs: Accounting for Periodic Boundary conditions

• Periodic cell-pair offsets translate particles once, avoiding per-pair wrapping conditions and additions in the inner loop.The stored offset is ±L for pairs requiring wrapping and zero otherwise.Source: §4.6 · p. 11

4.7 Speeding up the Calculation of Separations in Cell-Pairs: Sorting to enable Late Entry and Early Exit Conditions

• Sorting particles by z enables late entry and early exit in the j-loop, reducing the search volume from 3πmax to 2πmax.For each i-particle, the loop starts when dz > −Rmax and stops when dz ≥ Rmax.Source: §4.7 · p. 11

4.8 Speeding up the Calculation of Separations in Cell-Pairs: Explicit Vectorisation

• Explicit SIMD intrinsics vectorize j-particle separation calculations, and macros let one source file support float32 and float64 arrays.The macros expand at compile time into separate precision-specific functions, with the interface selecting the appropriate function for the input array type.Source: §4.8 · p. 12

4.9 Speeding up the Calculation of Separations in Cell-Pairs: Finding and Updating the Histogram Bins

• Squared distances and cos θ avoid sqrt and arccos for bin lookup, while reverse bin traversal makes runtime depend weakly on bin count for fixed Rmax.The reverse traversal reduces histogram-loop executions by checking later bins first.Source: §4.9 · pp. 12–13

4.10 OpenMP— Parallelising for Multi-cores

• Dynamic OpenMP scheduling assigns nearby primary cells to threads to improve cache reuse; the authors report very small load imbalance in practice.Each thread maintains a private histogram during computation, and these histograms are combined afterward; domain partitioning itself is not parallelized.Source: §4.10 · p. 13

4.11 Summary of optimisations in Corrfunc

• Corrfunc combines algorithmic reductions in distance computations with software optimizations that make each computation more efficient.The cited examples include partitioning, avoiding duplicate auto-correlation work, contiguous particle storage, explicit vectorization, and popcnt-based histogram updates.Source: §4.11 · p. 13

5.1 Complexity of the Code

• Cell-based searching gives O(NM) work, with typical Rmax values implying a 125–1000× speedup over brute force, but Corrfunc remains O(N2) in particle count at fixed Rmax.The code also uses custom approximations for division and arccos in DDrppi_mocks and DDtheta_mocks.Source: §5.1 · pp. 13–14

5.2 Scaling with Number of Particles

• For large particle counts, runtime scales as O(N2); at fewer than 10^4 particles, SIMD runtimes for wp(rp) and DR(θ) can decrease as particle count increases.The authors suspect more work shifts from the scalar remainder loop to efficient SIMD instructions as particle counts rise.Source: §5.2 · p. 14

5.3 Scaling with Maximum Search Radius

• Runtime recovers the expected O(Rmax3) and O(θmax2) scaling at larger search radii, after flat or decreasing runtimes at low Rmax.The authors speculate that low cell occupancy leaves more computation in the scalar remainder loop; larger occupancy shifts work to SIMD instructions.Source: §5.3 · p. 14

5.4 Scaling with OpenMP threads

• Auto-correlations reach 93% efficiency at 24 threads, while cross-correlations fall to 60–80% at 24 threads as cache utilization worsens.Cross-correlations remain 90% efficient to approximately 10 threads.Source: §5.4 · p. 14

5.5 Speedup from SIMD code

• AVX and SSE kernels deliver approximately 3× and 2× speedups, respectively, over the Fallback kernel.These gains are below the theoretical maxima, leaving scope for further vector-kernel optimization.Source: §5.5 · p. 14

6 COMPARISON WITH OTHER CODES

• The benchmark compares each code only with its corresponding Corrfunc runtime, checks pair-count agreement, and cautions that even unwitting setup choices may bias results.Tests were repeated three times and averaged; the authors encourage readers and code authors to benchmark independently.Source: §6 · p. 15
• In Rmax scaling tests, SSE and AVX kernels are progressively faster than the generic Fallback kernel; at 24 threads, efficiencies are 93%–105% for theory auto-correlations and 60%–80% for mock cross-correlations.Source: Figure 5 · p. 16, Figure 6 · p. 17

6.3 kdcount, version 0.3.21 (Feng et al. 2017)

• 3.1–6.0× single-threaded and 6.8–130.5× multi-threaded are the reported Corrfunc speed-ups over kdcount.The comparison used Corrfunc with autocorr = False.Source: §6.3 · p. 17

6.4 halotools, version 0.4 (Hearin et al. 2017)

• 1.3–4.8× single-threaded and 2.0–8.3× multi-threaded are the reported Corrfunc speed-ups over halotools.Fake RR counts were passed to avoid computing RR, and downsampling was disabled.Source: §6.4 · p. 17

6.5 TreeCorr, version 3.3.6 (Jarvis et al. 2004)

• 1.9–5.7× single-threaded and 1.1–7.3× multi-threaded are the reported Corrfunc speed-ups over TreeCorr.TreeCorr's bin_slop was disabled so its pair counts exactly agreed with Corrfunc.Source: §6.5 · p. 17

6.6 CUTE_box, git commit ab33dd8, (Alonso 2012)

• 1.6–6.5× single-threaded and 0.3–4.8× multi-threaded are the reported Corrfunc speed-ups over CUTE_box.CUTE_box was the only code with different pair counts; a small number of pairs appeared shifted by one bin, possibly because of floating-point effects at bin boundaries.Source: §6.6 · p. 17

6.7 mlpack RangeSearch, version 2.0.1 (Curtin et al. 2013)

• 0.8–8.7× is the reported Corrfunc speed-up over mlpack RangeSearch; this test reduced Rmax to 36 h−1Mpc because mlpack supports one bin and runs out of memory for many pairs.The reported mlpack timing sums tree-building and range-search times, which the authors say was significantly faster than its actual runtime, likely because it explicitly constructs and outputs every pair.Source: §6.7 · p. 18

6 COMPARISON WITH OTHER CODES

• Corrfunc is a factor of a few faster than other codes at moderate to high particle loads, but mlpack RangeSearch and CUTE_box outperform it below about 10^5 particles.At those low loads, the slower Corrfunc runtimes are fractions of a second and are described as unlikely to become a computational bottleneck.Source: §6.8 · p. 18

7 DISCUSSION

• Partitioning space into cells and searching neighboring cells reduces the search volume and the number of distance computations.Particles are sorted within cells to support short-circuit checks based on trigonometric inequalities.Source: §7 · p. 18
• Corrfunc beats the tested codes except mlpack RangeSearch at low particle numbers in the single-threaded benchmark, and except at very small particle numbers with 24 threads.The multi-threaded figure uses a different y-axis scale from the single-threaded figure.Source: Figure 7 · p. 18, Figure 8 · p. 18
• Corrfunc duplicates particle data for contiguous memory access, uses explicit SIMD instructions, batches histogram updates, and parallelizes pair counters with OpenMP.Batching groups pairs that fall in the same bin before incrementing that bin.Source: §7 · pp. 18–19
• Further proposed optimizations include refining cells and using 64-bit histogram processing for unweighted counts when average pair separation is not requested.The 64-bit approach is proposed to reduce calls to popcnt.Source: §7 · p. 19
• Loop-blocking may improve data reuse but initially slowed low-to-moderate particle cases; matrix-like kernels could also use BLAS or GPU libraries.The authors note that existing early-exit and late-entry checks may limit loop-blocking benefits.Source: §7 · p. 19
• Duplicated boilerplate across pair counters increases maintenance overhead and can be a source of bugs; autogenerated code is a proposed future direction.The common setup, cell-pair association, SIMD dispatch, and thread-result collection account for much of the duplicated code.Source: §7 · p. 19

8 CONCLUSIONS

• Corrfunc’s performance combines domain-based search pruning, cache locality, explicit vectorization, and OpenMP multicore parallelism.The authors report excellent strong scaling in their tests.Source: §8 · p. 19
• Spatial grids restrict distance calculations to pairs that could fall within Rmax, reducing the search volume.Particles are arranged contiguously within each cell to improve memory access speeds.Source: §8 · p. 19
• Explicit vector intrinsics compute multiple pairs simultaneously because the tested compilers did not generate vectorized correlation-function code.Source: §8 · p. 19
• Corrfunc uses OpenMP without shared mutable resources, and the authors report excellent strong scaling in their tests.Source: §8 · p. 19
• Designing, writing, and maintaining Corrfunc has been time-intensive, with its algorithms and source code evolving over more than five years.Source: §8 · p. 19

A1 Pair-counting in Cartesian volumes

• For Cartesian volumes, point separation uses the standard Euclidean distance.Source: A1 Pair-counting in Cartesian volumes · p. 21
• Corrfunc takes spatial-bin definitions from the first two columns of a text file and assumes z is the line of sight, with redshift-space distortions already applied to positions.The file is specified by filename on the command line or by a NumPy array through Python.Source: A1 Pair-counting in Cartesian volumes · p. 21
• DD(r) and ξ(r) use full three-dimensional separation r, with the first two bin-file columns specifying r.Source: A1 Pair-counting in Cartesian volumes · p. 21
• DD(rp, π) bins projected separation rp from the first two bin-file columns and uses linearly spaced π bins of width 1 up to πmax.The appendix gives π = |z1 − z2|.Source: A1 Pair-counting in Cartesian volumes · p. 21
• wp(rp) uses rp bins from the first two bin-file columns and includes all points with line-of-sight separation up to πmax.Source: A1 Pair-counting in Cartesian volumes · p. 21

• On-sky line-of-sight vector ℓ connects the observer to the midpoint of the segment joining the two points.Source: A2 Pair-counting in Spherical volumes · p. 21
• Figure A1 follows Fisher et al. (1994) conventions for defining separations between points on the sky.Source: Figure A1 · p. 21
• Projected on-sky correlation functions follow Fisher et al. (1994); rp bins occupy the first two file columns, while π bins have unit width through πmax.Source: A2 Pair-counting in Spherical volumes · p. 21
• Angular separation ∆σ is the angle between vectors v1 and v2, and Corrfunc handles angles from 0° to 180°.Source: A2 Pair-counting in Spherical volumes · p. 21

• When the first bin starts at zero, the auto-correlation DD term includes self-pairs by adding N to the appropriate bin after counting.Source: Appendix B · pp. 21–22
• The RR term uses particle density (N−1)/V because each primary particle has only N−1 other particles available for nonzero-distance pairs.Source: Appendix B · p. 22
• Using density N/V biases the estimator: a uniform random particle set yields ξ(r)=−1/N rather than 0.Source: Appendix B · p. 22
• For consistency with DD, RR must include self-pairs when the first bin includes zero separation.Source: Appendix B · p. 22