🔭 astrophysics

Fast Graph-based Higher-Order Clustering Statistics on the GPU

This paper presents a significant update to the GRAMSCI algorithm, introducing a GPU-accelerated port with optimized merge-walk sorting, parity-decomposed 4-point correlation functions, and out-of-core capabilities that achieve up to 9x speedups over CPU implementations for computing higher-order clustering statistics in cosmological datasets.

Original authors: Cristiano G. Sabiu

Published 2026-07-09

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Cristiano G. Sabiu

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to understand the shape of the universe by looking at how galaxies are scattered across the night sky. For a long time, astronomers have mostly looked at pairs of galaxies to see how often they appear close together. This is like looking at a crowd of people and only counting how many pairs of friends are standing next to each other. It tells you a lot, but it misses the bigger picture of how groups of three, four, or more people are interacting.

This paper introduces a major upgrade to a computer tool called GRAMSCI (which stands for "GRAph Made Statistics for Cosmological Information"). Think of GRAMSCI as a super-fast accountant that doesn't just count pairs, but counts every possible group of three, four, or more galaxies to reveal the hidden "clumping" patterns of the universe.

Here is what the authors did to make this tool faster and smarter, explained through everyday analogies:

1. The "Merge-Walk" Shortcut (Speeding Up the Count)

The Old Way: Imagine you have a list of neighbors for every person in a crowd. To find if Person A and Person B are both friends with Person C, the old computer had to look up Person C's list, search for Person A, then search for Person B. This was like looking up a name in a phone book over and over again. It worked, but it was slow.

The New Way: The authors realized that if you sort the lists alphabetically (or by ID number), you don't need to search. Instead, you can use two fingers to walk down both lists at the same time. If the names match, you count it! If one list is "ahead," you just move that finger forward.

The Result: This "merge-walk" method is much faster. For groups of three galaxies, it's about 1.5 to 2 times faster. For groups of four, it's nearly 4 times faster.

2. The "Left-Handed vs. Right-Handed" Detector (Parity Decomposition)

The Concept: In physics, some things look the same in a mirror, and some don't. A tetrahedron (a pyramid with a triangular base made of four galaxies) can be "left-handed" or "right-handed," just like your hands. You can't rotate a left hand to look like a right hand.

The Innovation: The new tool can now tell the difference between these "handed" groups. It separates the data into two buckets: "Even" (symmetric) and "Odd" (asymmetric).
Why it matters: If the universe is perfectly symmetric, the "Odd" bucket should be empty. If the tool finds a non-zero number in the "Odd" bucket, it suggests the universe might have a fundamental "handedness" or bias, which would be a huge discovery. The tool is smart enough to ignore "flat" shapes that don't have a clear handedness to avoid computer errors.

3. The "True Connection" Filter (Connected vs. Disconnected)

The Problem: When you count groups of four galaxies, many of those groups aren't actually "connected" in a special way. They are just two separate pairs of friends standing near each other by coincidence.

The Solution: The new tool automatically calculates what those "accidental" groups would look like and subtracts them.
The Analogy: Imagine you are counting how many people are in a "four-person huddle." The old tool counted everyone standing in a group of four, even if it was just two couples chatting. The new tool subtracts the two couples to find the actual four-person huddle. This reveals the "true" four-point signal that tells us about the early universe.

4. The "Super-Powered" GPU Engine

The Challenge: Counting these groups for millions of galaxies creates a massive amount of data. Doing this on a standard computer is like trying to move a mountain of sand with a teaspoon.

The Upgrade: The authors moved the entire calculation engine to a Graphics Processing Unit (GPU). Think of a CPU (standard computer brain) as a single master chef who is very smart but can only chop one vegetable at a time. A GPU is like a kitchen with 10,000 sous-chefs who can chop vegetables simultaneously.
The Result: The new tool is 2.6 to 9 times faster than the best 64-core computer.
The "Out-of-Core" Trick: Sometimes the data is so huge (like 45 gigabytes of galaxy connections) that it won't fit in the GPU's memory (which might only hold 24 GB). The authors built a "tiling" system. It's like moving a giant pile of bricks into a small room: you bring in a few bricks, do the work, move them out, and bring in the next batch. This allows the tool to handle data sets that are physically larger than the computer's memory.

5. Real-World Testing

The authors didn't just build this in theory; they tested it on real data from the DESI (Dark Energy Spectroscopic Instrument) survey, which maps millions of galaxies.

They compared their results against computer simulations (mocks) and found they matched perfectly.
They successfully measured the "Baryon Acoustic Oscillation" (a cosmic ruler left over from the Big Bang) using these complex group statistics, proving the tool works on real, messy, real-world data.

Summary

In short, the authors took a powerful tool for mapping the universe and made it faster (by changing how it searches lists), smarter (by detecting "handedness" and filtering out accidental groups), and stronger (by running on a super-fast graphics card that can handle data larger than its own memory). This allows scientists to ask deeper questions about the structure of the universe using data that was previously too difficult to process.

Technical Summary: Fast Graph-based Higher-Order Clustering Statistics on the GPU

Problem Statement
The spatial distribution of galaxies encodes critical information about the primordial density field and gravitational evolution. While the two-point correlation function (2pCF) has long been the standard tool for analyzing this distribution, it is insufficient for capturing the non-Gaussian structure generated by late-time nonlinear gravitational evolution and galaxy bias. Higher-order statistics, specifically the three-point (3pCF) and four-point (4pCF) correlation functions, provide complementary information, such as breaking degeneracies in bias parameters and probing primordial non-Gaussianity and parity violation. However, computing these statistics for large cosmological catalogues is computationally prohibitive. A brute-force approach scales as $O(N^k)$ , and while tree-based methods reduce this to $O(N m^{k-1})$ (where $m$ is the mean neighbor count), the cost remains high for $k=3$ and $k=4$ on modern survey scales (e.g., DESI). Furthermore, existing codes often lack direct configuration-space implementations of connected 4pCFs with parity decomposition or efficient GPU acceleration capable of handling graphs exceeding device memory.

Methodology
This paper presents Version 2 of GRAMSCI (GRAph Made Statistics for Cosmological Information), an algorithm that utilizes a graph-database approach built upon a $k$ -d tree. The core methodology involves constructing a compact, sorted adjacency structure from the initial coordinates, which is then queried to count $N$ -tuples without revisiting the raw data.

The update introduces five key methodological advancements:

Merge-Walk Algorithm: The original binary-search inner loop, which incurred an $O(m \log m)$ cost per hub-spoke pair, is replaced by a merge-walk algorithm. By exploiting the sorted order of adjacency lists, the algorithm performs a linear pass to find set intersections, reducing the complexity to $O(m)$ .
Parity-Decomposed 4pCF: A new estimator separates the 4pCF signal into parity-even ( $\zeta_+$ ) and parity-odd ( $\zeta_-$ ) channels. This is achieved by calculating the signed volume of galaxy tetrahedra. To avoid floating-point noise contamination in the parity-odd channel, tetrahedra with degenerate volumes (effectively zero) are deterministically assigned zero weight in the odd channel.
Connected 4pCF Estimation: The code now internally computes the disconnected contribution (products of 2pCFs) on the same graph and binning. This allows for the direct subtraction of the disconnected term from the total 4pCF to isolate the connected (trispectrum) signal, which is essential for studying primordial non-Gaussianity.
Python Interface: A lightweight wrapper allows the Fortran engine to be called directly from NumPy arrays, handling file I/O and returning structured result objects.
GPU Port with Out-of-Core Support: The complete query engine (3pCF, 4pCF, and parity 4pCF) has been ported to NVIDIA GPUs using OpenACC. The implementation uses a Compressed-Sparse-Row (CSR) data layout and employs a "gang-and-vector" mapping strategy where one thread block processes one hub, and vector lanes handle the inner neighbor loops. To handle graphs larger than device memory (e.g., BAO-scale surveys), an out-of-core tiling scheme splits the edge arrays into windows, processing the graph in chunks while maintaining exact counting.

Key Results
The authors benchmarked the new code against the original CPU implementation and analytic injections:

Algorithmic Speedup: The merge-walk kernel provides a speedup of $1.5\times$ to $2.0\times$ for the 3pCF and $2.8\times$ to $3.8\times$ for the 4pCF over the binary-search approach on a 64-thread CPU node, with gains increasing as the maximum scale ( $r_{max}$ ) and neighbor count grow.
GPU Performance: On a single consumer GPU (NVIDIA RTX 3090 Ti), the code achieves speedups of $2.6\times$ (3pCF) to $9\times$ (4pCF) compared to a 64-thread CPU node. The GPU performance scales with memory bandwidth rather than floating-point throughput.
Scalability: The out-of-core tiling scheme allows the processing of graphs far exceeding device memory. A specific test measured a $9 \times 10^9$ -edge 3pCF on a 24 GB card with only $\sim20\%$ overhead compared to a hypothetical all-resident run.
Correctness: The GPU kernels reproduce the CPU reference results exactly in binned counts. Validation against analytic signal-injection tests for parity violation and null tests on parity-symmetric simulations (Quijote) confirms the estimator is unbiased and correctly handles degenerate geometries.
Application: The code was applied to the DESI DR1 LRG sample (3.0 million galaxies). The measured 3pCF and connected 4pCF track the EZmock ensemble predictions, with the 3pCF showing the expected BAO feature and the connected 4pCF isolating the non-Gaussian signal.

Significance
The paper claims that these updates significantly expand the practical utility of higher-order statistics in cosmology. By reducing the computational cost of the 4pCF by nearly an order of magnitude on consumer hardware, the code makes it feasible to evaluate connected 4pCFs across the full BAO range for large catalogues like DESI, a task previously limited by the need for thousands of mock realizations for covariance estimation. The inclusion of a parity-decomposed estimator provides a direct, configuration-space cross-check for parity violation searches, independent of spherical-harmonic analyses. Finally, the ability to compute connected 4pCFs internally without external inputs ensures self-consistency in binning and pair sets, which is critical for isolating the non-Gaussian trispectrum signal. The authors conclude that the code is publicly available and ready for application to current and future large-scale structure surveys.