Galactos: Computing the Anisotropic 3-Point Correlation Function for 2 Billion Galaxies

The nature of dark energy and the complete theory of gravity are two central questions currently facing cosmology. A vital tool for addressing them is the 3-point correlation function (3PCF), which probes deviations from a spatially random distribution of galaxies. However, the 3PCF’s formidable computational expense has prevented its application to astronomical surveys comprising millions to billions of galaxies. We present Galactos, a high-performance implementation of a novel, O(N^2) algorithm that uses a load-balanced k-d tree and spherical harmonic expansions to compute the anisotropic 3PCF. Our implementation is optimized for the Intel Xeon Phi architecture, exploiting SIMD parallelism, instruction and thread concurrency, and significant L1 and L2 cache reuse, reaching 39% of peak performance on a single node. Galactos scales to the full Cori system, achieving 9.8PF (peak) and 5.06PF (sustained) across 9636 nodes, making the 3PCF easily computable for all galaxies in the observable universe.

💡 Research Summary

The paper introduces Galactos, a high‑performance implementation designed to compute the anisotropic three‑point correlation function (3PCF) for galaxy surveys containing up to two billion objects. The 3PCF encodes richer cosmological information than the two‑point correlation function (2PCF), including sensitivity to galaxy bias, non‑linear growth, and redshift‑space distortions, which are crucial for probing dark energy and testing modifications to General Relativity. Historically, the 3PCF has been computationally prohibitive because a naïve algorithm scales as O(N³), making it infeasible for modern surveys.

Galactos overcomes this barrier by employing a novel O(N²) algorithm that builds on the Legendre‑polynomial approach used for the isotropic 3PCF but extends it to capture anisotropies. The method treats each galaxy as a “primary” and gathers all secondary galaxies within a maximum radius (Rmax = 200 Mpc/h) using a load‑balanced k‑d tree. For each radial shell around the primary, the density field is expanded in spherical harmonics Yℓm. By exploiting the spherical‑harmonic addition theorem, the algorithm combines the coefficients from two shells to obtain the multipole moments ζℓm(r₁,r₂) that fully describe the anisotropic 3PCF. Because the expansion decouples angular dependence from radial binning, the dominant computation reduces to pairwise operations, yielding O(N²) complexity.

Implementation details are tuned for the Intel Xeon Phi (Knights Landing) architecture. The authors reorganize data into a structure‑of‑arrays layout to enable contiguous memory accesses, and they vectorize the core kernels with 512‑bit AVX‑512 instructions, achieving about 39 % of the theoretical peak FLOP rate per node. Thread‑level parallelism is managed with OpenMP, while MPI distributes the k‑d tree partitions across 9 636 nodes of the Cori supercomputer. The design ensures near‑perfect weak and strong scaling, with a reported peak performance of 9.8 PF and sustained performance of 5.06 PF.

Performance benchmarks demonstrate that Galactos can compute the full anisotropic 3PCF for a realistic 2‑billion‑galaxy catalog in roughly two hours, a task that would be impossible with previous O(N³) methods. Compared to state‑of‑the‑art 2PCF implementations, Galactos delivers a five‑fold speedup on the same hardware while providing substantially more cosmological information.

The authors conclude that the combination of algorithmic innovation (O(N²) anisotropic expansion), aggressive SIMD/vector optimization, and scalable k‑d tree load balancing makes the 3PCF a practical tool for upcoming large‑scale redshift surveys such as DESI, Euclid, and LSST. By enabling routine measurement of the anisotropic 3PCF, Galactos opens a new avenue for tightening constraints on dark energy, testing gravity theories, and improving our understanding of galaxy formation.