arXiv CSopen access

Preserving Clusters in Error-Bounded Lossy Compression of Particle Data

Preserving Clusters in Error-Bounded Lossy Compression of Particle Data

arXiv:2604.18801v1 [cs.LG] 20 Apr 2026

Congrong Ren The Ohio State University Columbus, Ohio [email protected]

Sheng Di Argonne National Laboratory Lemont, Illinois [email protected]

Franck Cappello Argonne National Laboratory Lemont, Illinois [email protected]

Abstract—Lossy compression is widely used to reduce storage and I/O costs for large-scale particle datasets in scientific applications such as cosmology, molecular dynamics, and fluid dynamics, where clustering structures (e.g., single-linkage or Friends-ofFriends) are critical for downstream analysis; however, existing compressors typically provide only pointwise error bounds on particle positions and offer no guarantees on preserving clustering outcomes, and even small perturbations can alter cluster connectivity and compromise scientific validity. We propose a correction-based technique to preserve single-linkage clustering under lossy compression, operating on decompressed data from off-the-shelf compressors such as SZ3 and Draco. Our key contributions are threefold: (1) a clustering-aware correction algorithm that identifies vulnerable particle pairs via spatial partitioning and local neighborhood search; (2) an optimizationbased formulation that enforces clustering consistency using projected gradient descent with a loss that encodes pairwise distance violations; and (3) a scalable GPU-accelerated and distributed implementation for large-scale datasets. Experiments on cosmology and molecular dynamics datasets show that our method effectively preserves clustering results while maintaining competitive compression performance compared with SZ3, ZFP, Draco, LCP, and space-filling-curve-based schemes. Index Terms—lossy compression, error control, scientific data, clustering

I. I NTRODUCTION Large-scale particle simulations represent a dominant workload in modern high-performance computing (HPC), where error-bounded lossy compression has emerged as a vital data reduction strategy to mitigate the immense storage, transfer, and analysis overhead. In cosmology, the Hardware/Hybrid Accelerated Cosmology Code (HACC) framework [1] is currently scaling toward tens of trillions of particles to model the universe’s evolution. At the Exascale, these simulations generate individual snapshots exceeding 500 TB, with aggregate I/O throughput peaking at over 30 TB/s and cumulative data products reaching the exabyte scale [2]. Similarly, molecular dynamics simulations in biology and materials science now produce trillion-particle snapshots [3] to study complex phenomena such as polymer clustering [4] and shock-induced plasticity [5]. Unlike lossless compression, which preserves the

Katrin Heitmann Argonne National Laboratory Lemont, Illinois [email protected]

Hanqi Guo The Ohio State University Columbus, Ohio [email protected]

data bit-for-bit but yields limited compression ratios (around 2:1 for floating-point data [6]), error-bounded lossy compression significantly reduces the data volume while guaranteeing high data quality, and thus has been widely adopted across scientific domains including cosmology [7], [8], materials science [9], and fluid dynamics [10], [11]. Despite their success, existing lossy compression methods face a significant limitation when applied to particle data: they only guarantee pointwise error bounds and may inadvertently destroy cluster memberships, leading to distorted structural relationships and incorrect scientific inferences. General-purpose scientific compressors, such as SZ3 [12] and ZFP [13], are primarily designed for regular Cartesian grids; when applied to discrete, irregular particle distributions, they often yield suboptimal compression ratios because they do not exploit the spatial coherency of particles. While more recent compressors like Google Draco [14] and LCP [9] are specifically tailored for particle data, they still focus exclusively on bounding coordinate reconstruction errors. For many scientific workflows, however, bounding pointwise error is insufficient for maintaining the integrity of downstream structural analyses. A prominent example is single-linkage clustering, which underlies the friends-of-friends halo finder in cosmology [15], [16], the detection of coherent structures in fluid dynamics [17]– [19], and the identification of ion clusters [20] and protein aggregation [21] in molecular simulations, where the precise determination of whether particle pairs lie within a given distance threshold is critical. Currently, cluster invariance can only be ensured indirectly. One approach is to impose strict enough global error bounds to avoid artificial merging or splitting of clusters. However, this often necessitates near-lossless thresholds, which severely limits achievable compression ratios. Our empirical tests on an HACC dataset [22] with over one billion 3D particles demonstrate that SZ3 and LCP require a relative error bound around 10−8 to preserve cluster membership and attain compression ratios of only 1.42 and 1.75, respectively. These low ratios arise from treating all particles as equally important, ignoring

that most coordinates do not affect clustering outcomes and could be stored at lower precision. An alternative approach, storing cluster assignments as a separate integer array alongside the lossy-compressed coordinates, decouples cluster preservation from global error bounds but incurs significant auxiliary storage overhead. While cosmological analyses often ignore the clusters below a certain particle threshold, the indices and cluster assignment of approximately 0.1N −0.2N of particles still require storage, where N is the number of particles. For example, the HACC dataset [22] with 1.07 billion particles has 113.8 and 168.7 million particles within clusters meeting typical size thresholds of 100 [23], [24] and 20 [25], respectively; storing these assignments using 4-byte integers requires an additional 7.1% and 10.5% of its original storage. Furthermore, this overhead multiplies when users explore the parameter space of clustering, as each parameter setting (e.g., linking length) requires a unique assignment array. Also, this metadata-heavy approach does not improve the fidelity of the compressed coordinates themselves, leaving other downstream analyses unimproved. These limitations necessitate methods that achieve high cluster fidelity while maintaining good compression ratios. We propose a novel correction-based method that preserves single-linkage cluster membership under lossy compression with a user-defined error bound ξ and distance threshold b. The method is applied in situ during compression and is compatible with any off-the-shelf lossy compressor (referred to as the base compressor). It efficiently detects and corrects broken particle pairs near the clustering threshold, which are identified using a spatial cell-based partitioning technique. Specifically, we formulate cluster preservation as a constrained optimization problem, where the loss function penalizes violations of pairwise distance constraints. We then apply projected gradient descent to restore the original connectivity graph while strictly preserving the error bound enforced by the base compressor. To ensure scalability, we implement a data-parallel solver with both CPU and GPU backends, supporting single-node systems and multi-node MPI environments. We evaluate our solution using large-scale simulation datasets across multiple application domains including cosmology, molecular dynamics, and fluid dynamics. Compared to state-of-the-art error-bounded compressors (SZ3, ZFP) and geometry-aware methods (Draco, LCP), our correction technique consistently achieves superior preservation of singlelinkage clusters with negligible storage overhead and computational cost. We evaluate clustering accuracy using the Matthews correlation coefficient (MCC) for pairwise link existence and the halo mass function (HMF) for particle cluster assignments. Furthermore, performance analysis demonstrates that our parallel GPU solver achieves up to 62× speedup over the CPU baseline, confirming that the correction step integrates into high-throughput simulation pipelines without introducing a computational bottleneck. Our contributions are three-fold: • A Novel Correction Algorithm: We propose a method for preserving the single-linkage clusters of particle data while satisfying user-defined coordinate error bounds;

GPU and Distributed Parallelism: We design and implement a distributed, data-parallel solver that scales across multi-node CPU and GPU clusters via MPI; • Comprehensive Evaluation: We provide a rigorous validation across multiple scientific datasets and state-of-theart compressors using clustering metrics MCC and HMF. •

II. R ELATED W ORK AND BACKGROUND We summarize literature on error-bounded lossy compression, especially feature-preserving compression, alongside scientific applications of single-linkage clustering. A. Error-bounded Lossy Compression Error-bounded lossy compressors ensure that the pointwise difference between the original and decompressed values remains within a user-specified bound. Lossy compression for regular grid data is broadly categorized by their reconstruction primitives [26]. For example, prediction-based methods, such as SZ family including SZ3 [12] and cuSZp2 [27], use predictors (e.g., the Lorenzo predictor [28]) to estimate data values based on their neighbors, subsequently quantizing the residual error to satisfy strict L∞ guarantees. In contrast, transform-based methods, such as ZFP [13], employ block-based transforms to map data into sparse coefficients that are more amenable to compression. While these tools achieve high compression ratios for smooth, structured fields, they rely on the implicit spatial locality of grid indices. For irregular data, where index-level adjacency does not correspond to spatial proximity, these methods yield significantly degraded compression ratios. Lossy compression for irregular data, such as unstructured meshes and gridless particle data, uses spatial distribution or proximity to restore the locality lost in flat indexing. For example, Google Draco [14] employs octree-based partitioning and coordinate quantization for efficient compression, while the MPEG G-PCC [29] and V-PCC [30] standards similarly apply octree-based spatial decomposition to target geometryand video-based point clouds. Domain-specific solutions like MDZ [31] exploit unique spatial patterns (e.g., zigzag or stairwise) and high temporal similarity in molecular dynamics trajectories. More recent advancements, such as LCP [9], use block-wise decomposition to improve data-parallel processing. However, these methods remain fundamentally “featureagnostic:” they prioritize pointwise accuracy but do not explicitly guarantee the preservation of derived clusters. Feature-preserving lossy compression addresses the limitations of pointwise error bounding by ensuring that specific quantities of interest (QoIs) remain accurate after decompression. These QoIs typically include structural, statistical, physical, and topological properties that are critical for downstream scientific tasks and analyses [26]. For example, QPET [32] is designed to maintain the quantile distribution of the data. Post-compression correction-based (also referred to as editor augmentation-based) approaches have emerged as an effective strategy to satisfy QoI requirements. Instead of modifying the base compression algorithm itself, these methods refine the

decompressed output by applying “edits” to restore specific features while strictly adhering to the original pointwise error bounds. For example, to preserve power spectra, FFCz [8] derives edits by alternating projections between spatial and frequency constraints. To preserve scalar field topology, MSz [33] edits decompressed data to maintain Morse–Smale segmentations, while Gorski et al. [34] focus on contour tree accuracy. However, unlike these existing methods which focus on regular-grid scalar fields, preserving the integrity of singlelinkage clusters in particle data remains an open challenge, as it requires managing the high sensitivity of distance to small coordinate changes near the linking threshold. B. Single-linkage clustering Single-linkage clustering, known as the friends-of-friends (FoF) algorithm in cosmology, is a foundational tool across numerous scientific domains [15], [17], [20]. This method defines clusters based on a proximity rule: for a set of N −1 particles P = {pn }N n=0 , the FoF algorithm identifies clusters based on a proximity threshold called the linking length denoted by b. We define a connectivity graph G = (P, E) where an edge (pi , pj ) ∈ E exists if and only if the Euclidean distance of the two particles pi and pj satisfies d(pi , pj ) ≤ b. Each FoF cluster (referred to as a “halo” in cosmology [35]) corresponds to a connected component of G (Fig. 1).

(a)

(b)

(c)

Fig. 1. Illustration of the FoF clustering algorithm. (a) An initial distribution of particles in space. (b) Edge construction for all particle pairs with a Euclidean distance d(pi , pj ) ≤ b, where b is the linking length. (c) The resulting connected components, which define the final clusters (halos).

In practice, b is determined by a linking parameter η and the mean particle separation ∆p . For a dataset in d-dimensional space with a total volume Vvol , b is defined as:  1/d Vvol b = η∆p , where ∆p = N In cosmological simulations, η is typically set to 0.2 or 0.168 [1], [15], [36]. This formulation ensures that the clustering threshold scales appropriately with the global density. 1) Single-linkage Clustering Implementation: Efficiently computing FoF clusters for large-scale datasets requires specialized spatial indexing and parallel algorithms to avoid the O(N 2 ) complexity of a naive pairwise search [1], [37]. A widely adopted strategy involves spatial partitioning and cell linking, where the simulation volume is partitioned into a uniform grid of cells [1]. By setting the cell side length w ≥ b, the search for any particle pair (pi , pj ) such that d(pi , pj ) ≤ b is restricted to pairs located within the same cell or two adjacent

cells. In parallel environments, each thread typically processes a unique cell and evaluates candidate pairs within that cell and its forward neighboring cells. By only checking neighbors in one direction (e.g., 4 specific neighbors in 2D or 13 in 3D), the algorithm avoids redundant distance computations for the same pair of cells. This cell-linking strategy reduces the complexity of edge identification to approximately O(N · ρ̄), where ρ̄ represents the average number of particles per cell. 2) Single-linkage Clustering Applications: In cosmology, FoF is a standard for halo finding, used to identify dark matter structures in N-body simulations [15], [16]. Similarly, in molecular dynamics, it is used to identify ion-pair clusters and characterize molecular self-assembly processes [20], [21]. It is also critical in fluid dynamics for detecting coherent structures in turbulent flows [17]–[19]. In these contexts, the resulting clusters, including mass functions and spatial distributions, form the essential basis for high-level scientific inference. The halo mass function (HMF), denoted as dn/d log M , represents the number density of dark matter halos per unit mass interval and serves as a critical downstream metric for validating cosmological simulations [38]. The HMF is a fundamental probe of the Λ Cold Dark Matter (ΛCDM) model, as its shape and evolution are directly governed by the growth of cosmic structures and the underlying dark matter density field [23], [39], [40]. For each identified halo i, we first calculate its total mass Mi ; assuming uniform particle mass, Mi ∝ Ni , where Ni is the number of constituent particles in halo i. We then partition these masses into B equal-width logarithmic bins and normalize the resulting counts by the total simulation volume and the logarithmic bin width. Maintaining high numerical precision in these calculations is vital, as FoF halo algorithm exhibits extreme sensitivity to coordinate perturbations near the linking threshold. Specifically, the addition of a single artificial link can merge two massive clusters, while the removal of a valid link can fragment a single halo into multiple components. Such inaccuracies in the HMF can lead to a fundamental mischaracterization of the universe, such as underestimating the total amount of matter or miscalculating how “clumpy” the cosmic structure is. III. M ETHODOLOGY This section formulates the problem of preserving FoF clusters as a constrained optimization problem, describes our iterative refinement algorithm, and details its implementation for single-node GPU and multi-node environments. While our current implementation assumes a non-periodic domain, the method extends straightforwardly to periodic boundary conditions by replacing Euclidean distances with minimum image distances in both vulnerable pair detection and loss computation, as is standard in cosmological N-body codes. A. Constrained optimization problem We formulate the preservation of FoF clusters as a constrained optimization problem, where the objective is to minimize the cumulative distance overflow of violated pairs subject

to the global pointwise error bound. Let the original 3D parti−1 cle data be P = {pn }N n=0 with coordinates (xn , yn , zn ), and the corresponding decompressed data from a base compressor −1 be P̂ = {p̂n }N n=0 with coordinates (x̂n , ŷn , ẑn ). Given an absolute error bound ξ, the reconstructed coordinates satisfy |x̂n − xn | ≤ ξ, |ŷn − yn | ≤ ξ, and |ẑn − z√n | ≤ ξ for all n, implying a maximum displacement of 3ξ for any single particle. Consequently, the Euclidean distance d(p̂i , p̂j ) between two arbitrary decompressed particles can deviate from √ √ its original value d(pi , pj ) by at most 2 3ξ. This 2 3ξ bound represents the worst-case scenario where two particles are displaced in opposite directions along the cube diagonals. While the average distance deviation is significantly lower in practice, this conservative radius is necessary to identify potential connectivity changes. Under a linking length √ b, links between pairs with an original distance in (b − 2 3ξ, b] may be inadvertently broken, while pairs whose distance is in √ (b, b + 2 3ξ] may result in falsely created links. √ √ We define pairs within the original distance range (b − 2 3ξ, b + 2 3ξ] as vulnerable pairs, as their connectivity is susceptible to compression errors. Accordingly, any particle belonging to at least one vulnerable pair is designated an editable particle, as modifying only their positions in the decompressed data is sufficient to restore the correct clusters. To preserve the clusters, we define a loss function L that penalizes distance deviations for pairs whose connectivity has been altered. The loss is formulated as: X X L(P̂ ) = (dˆi,j −b)2 + (b−dˆi,j )2 , (1) i,j:di,j ≤b<dˆi,j

i,j:dˆi,j ≤b<di,j

where di,j = d(pi , pj ) and dˆi,j = d(p̂i , p̂j ). The first term accounts for inadvertently broken links, while the second term accounts for falsely created links. By definition, only vulnerable pairs contribute to this loss, and only editable particles are involved in the optimization. Let E denote the set of editable particles. We then formulate the following constrained optimization problem to refine the particle positions: min

L(P̂ )

s.t.

|x̂n − xn | ≤ ξ,

p̂n ∈E

(2) |ŷn − yn | ≤ ξ,

|ẑn − zn | ≤ ξ

B. Projected gradient descent We propose a projected gradient descent (PGD) approach to iteratively refine particle coordinates, where each iteration performs a gradient descent to minimize the violation loss, followed by a projection of the updated coordinates back onto the box constraints defined by ξ. Alg. 1 shows the pseudocode. Identifying vulnerable pairs and editable particles. For a given linking length b and error bound ξ, we identify vulnerable pairs using a spatial partitioning and cell-linking strategy, as detailed in Section √ II-B. Specifically, we set the cell side length to w ≥ b + 2 3ξ, representing the maximum possible distance for a vulnerable pair. To maintain computational efficiency and avoid an overabundance of empty cells in sparse regions, we constrain the total number of cells to be no greater

Algorithm 1 FoF Cluster Preserving Correction Require: Original positions P , base-compressed positions P̂ (0) , error bound ξ, linking length b, bit depth m, learning rate α, max iterations Tmax , convergence threshold ϵL Ensure: Compressed bitstream of FoF cluster-preserved edits ϵq ← 2ξ/(2m − 1) ξ ′ ← ξ(1 − 2−m ) √ √ V ← {(i, j) : dij ∈ (b − 2 3ξ, b + 2 3ξ]} (0) P̂ ← P̂ for t = 1 to Tmax do if Ltight (P̂ ) ≤ ϵL then break end if g ← ∇P̂ Ltight (P̂ ) P̂ ← P̂ − α g P̂ ← projB(ξ′ ) (P̂ ) end for ∆ ← P̂ − P̂ (0) flags ← bitmask of non-zero entries in ∆, packed into 8-bit integers edits ← non-zero values of ∆, quantized to m bits return Huffman+ZSTD(flags, edits)

than the number of particles N . We then iterate through all particle pairs within the same cell or adjacent cells, using a directional search to avoid redundant cell pairs, to identify the set of vulnerable pairs and their constituent editable particles. Gradient descent and projection. The optimization begins by calculating the initial loss and terminating if it falls below a predefined tolerance. Otherwise, we enter an iterative refinement process. In each iteration, we first update the coordinates of the editable particles by gradient descent, and then apply a box projection to the updated coordinates. This ensures that every particle p̂n remains within its allowed error bound ξ, effectively enforcing the constraints defined in Equation (2). The process repeats until the loss satisfies the convergence criteria or the maximum iteration count is reached. Compaction, quantization, and lossless compression. Upon termination of the optimization, the coordinate shifts are recorded as “edits”, and decomposed into two distinct components: flags and compact edits. The flags consist of a bit-mask of length 3N that indicates which coordinates were modified; these are packed into 8-bit integers for efficient storage. The compact edits store only the non-zero values, the count of which is at most 3|E|. These edits undergo uniform quantization into 2m intervals, where m denotes the bit depth. To ensure that the quantized edits strictly satisfy the box constraints in Eq. (2), we apply a safety margin by shrinking the error bound to ξ ′ = ξ(1 − 2−m ). Furthermore, to ensure the restored clusters are robust against the maximum possible quantization error ϵq = 2m2ξ−1 , in practice we optimize a tightened loss Ltight in place of L in Eq. (2) during the PGD: X √ Ltight (P̂ ) = (dˆi,j − b + 2 3ϵq )2 (3) √ i,j:di,j ≤b,dˆi,j >b−2 3ϵq

+

X √ i,j:di,j >b,dˆi,j ≤b+2 3ϵq

√ (b + 2 3ϵq − dˆi,j )2 .

This refined loss over-corrects particle positions to create a safety zone that prevents bit-depth reduction from pushing distances back across the threshold b. Any pair corrected to

satisfy Ltight = 0 remains correctly classified after quantization, since the safety margin ϵq is chosen to exceed the maximum quantization error introduced by the edit encoding. To assess sensitivity to m, we evaluate m ∈ {8, 16, 32} on a molecular dynamics dataset, EXAALT, as a representative case; zero violations are re-introduced by quantization in all cases, confirming robustness across a wide range of bit-depths. Among these, m = 16 achieves the best compression ratio due to the tradeoff between edit precision and encoding overhead, and is used as the default in all experiments. Finally, the flags and quantized edits are compressed using Huffman coding [41] followed by ZSTD [42] to further reduce the storage. Reconstruction of the edited decompressed data. To reconstruct the data, we first apply ZSTD decompression and Huffman decoding to retrieve the flags and quantized edits. These compact edits are then dequantized and mapped back to their original coordinate indices using the bit-mask flags to reconstruct a full-length edit vector. The final coordinates are obtained by performing an element-wise addition of these reconstructed edits to the initial output of the base compressor. Convergence analysis. In practice the algorithm projects onto the tighter box B(ξ ′ ) with ξ ′ = ξ(1 − 2−m ) < ξ to absorb quantization error; the original positions P lie at ′ the center of B(ξ) and therefore also within √ B(ξ ), so P is ′ feasible within B(ξ ); moreover, since 2 3ϵq ≪ ξ, nearby positions achieve Ltight = 0, confirming the global minimum is zero. Although the loss is non-convex due to the falselink term, both components are C 1 with Lipschitz-continuous gradients since each term has the form [max(0, g)]2 where g is smooth [43], enabling monotonic loss decrease and O(1/T ) convergence to a stationary point under PGD. Two structural properties promote convergence to the global minimum in practice. First, false links among vulnerable pairs are rare, especially at small error bounds, so the loss is dominated by its convex broken-link component in most configurations. Second, the vulnerable-pair graph decomposes into many small independent components optimized separately; within each component, the geometric conditions required for gradient cancellation between conflicting pairs are extremely unlikely to be satisfied simultaneously. Together, these properties explain why PGD reaches L = 0 in all tested configurations when run to convergence without iteration cap, as confirmed empirically in Section IV, and the iteration budget scales with the number of violations rather than total particle count. The same analysis applies to Ltight , whose global √minimum of zero is achievable since the additional margin 2 3ϵq ≪ ξ is easily accommodated within the box constraints. In practice, we use Adam [44] for faster convergence; the theoretical guarantees above apply to vanilla PGD, and empirical convergence is confirmed in Section IV. C. Single GPU Acceleration Our GPU implementation parallelizes all stages. During spatial partitioning, each thread maps one particle to its grid cell and atomically increments per-cell counts in global memory; a device-wide exclusive prefix sum (via the CUB library [45])

then converts these counts into cell-start offsets, after which a scatter kernel places each particle into a cell-sorted array using atomic index reservations. This cell-sorted layout ensures that subsequent neighbor searches access spatially co-located particles, improving global memory coalescing. For vulnerable pair detection, the discovered pairs are buffered in thread-local register arrays and flushed to global memory in batches via a single atomic addition, significantly reducing contention on the global pair counter. Editable particles are registered in a direct-address map of size N in global memory; an atomicCAS on the particle’s index slot prevents duplicate insertions, yielding O(1) lookup during the PGD. In the PGD phase, threads compute per-pair loss and gradients, which are summed through a shared-memory reduction within each block and accumulated into a global buffer via a single atomic addition per block. Gradients are computed in a per-pair fashion and accumulated into per-coordinate buffers through global-memory atomic additions indexed by the direct-address map. The gradient update and box projection are then applied embarrassingly in parallel, with one thread per editable particle. D. Inter-node distributed-Memory Scaling To scale the algorithm to cosmological volumes exceeding a single GPU’s memory, we provide a distributed-memory implementation using the Message Passing Interface (MPI). The communication overhead consists of three phases. First, bounding box exchange is performed by first padding each √ rank’s bounding box by a ghost zone of width δ = b + 2 3ξ (in 3D) and then executing an MPI_Allgather to exchange bounding box metadata; two ranks are neighbors if their padded volumes overlap. Second, ghost exchange transfers particle coordinates from neighboring ranks via non-blocking MPI_Isend and MPI_Irecv. To hide latency, we overlap MPI transfers with GPU allocation of the extended particle buffer and copying local particles into it; execution pauses only for a synchronization wait before spatial partitioning begins. Third, global convergence check is monitored via an MPI_Allreduce of the loss function at every iteration, synchronizing all ranks before the subsequent gradient step. E. Complexity and scalability analysis Using a linked-list cell grid of cell side length w, each particle checks its own and 13 neighboring cells to identify vulnerable pairs, costing O(N ρ) where ρ is the average number of particles within the search region, identical to the cost of FoF itself. Each optimization iteration processes all |V | vulnerable pairs for gradient computation, giving O(T · |V |) total optimization cost. The iteration budget T ≤ ⌈12ξ 2 |Vtight |/ϵL ⌉ follows from bounding the initial loss Ltight (P̂ (0) ) by 12ξ 2 |Vtight |, where Vtight ⊇ Vviolated is the active pair set √ of Ltight ; each currently-violated pair contributes at most (2 3ξ)2 , and nearthreshold pairs not in Vviolated contribute at most 48ϵ2q ≪ 12ξ 2 , so |Vtight | ≈ |Vviolated | in practice. On the GPU, one thread is assigned per vulnerable pair and gradient accumulation uses a parallel tree reduction, reducing per-step complexity

to O(T log N ). In the distributed setting, each of the R MPI ranks communicates only boundary particles, with overhead O((N/R)2/3 ) per rank reflecting the surface-to-volume ratio of each spatial subdomain. IV. E VALUATION We present the evaluation scheme and summarize key findings in this section. A. Evaluation scheme Baselines and base compressors. We evaluate our method using SZ3 [12], ZFP [13], cuSZp2 [27], Google Draco [14], and LCP [9] as representative baselines and base compressors. SZ3 and ZFP are state-of-the-art general-purpose lossy compressors that provide strict pointwise error control and are widely used across diverse scientific applications. cuSZp2 is a GPU-accelerated compression technique based on SZ3. Draco and LCP are specifically designed for particle data. Furthermore, we implement a Morton-order baseline that operates post-spatial partitioning; this method traverses particles within each cell according to their Morton codes, applies a firstorder predictor and subsequent error quantization to provide a reference for spatial-locality-based compression. Metrics. We compare these baselines against the same base compressors corrected by our method using quantitative metrics, including compression ratio, rate distortion, and throughput. To evaluate the accuracy of the FoF clusters, we employ two specialized statistical metrics: the Matthews correlation coefficient (MCC), which measures the consistency of binary existence of connectivity of vulnerable pairs in the original and decompressed or edited data. It is defined as:

MCC = √

T P ×T N −F P ×F N , (T P +F P )(T P +F N )(T N +F P )(T N +F N )

where T P and T N count correctly linked and unlinked pairs, respectively; a false positive (F P ) occurs if a pair is unlinked in original data but becomes linked in decompressed data, while a false negative (F N ) represents a link in original data is broken in the decompressed data. An MCC of 1 denotes a perfect clustering match, while 0 indicates a correlation no better than random chance. We restrict MCC to vulnerable pairs, as non-vulnerable pairs are guaranteed to remain connected and would otherwise inflate the TN count, masking the true correction quality. Furthermore, to assess the impact on downstream cosmological analysis, we compare the HMFs derived from the original and reconstructed datasets. Datasets. The datasets used in our evaluation, summarized in Table I, span multiple scientific domains, spatial and temporal resolutions, and data modalities, and represents the diverse challenges faced by scientific data analysis. HACC data is generated from cosmological simulations carried out with the HACC framework [1], [46]. EXAALT data comes from molecular dynamics simulations [47]. Finite pointset method (FPM) datasets simulate the chemical process of salt dissolution in water, where higher-density salt diffuses into water to form unstable, downward-reaching structures known as viscous fingers [48]. For HACC (hiRes), results are reported

on a single representative timestep; for FPM datasets, all available timesteps are evaluated and results are averaged, as the method’s behavior is consistent across timesteps. TABLE I B ENCHMARK DATASETS . A LL DATASETS ARE IN 3D, SINGLE - PRECISION FLOATING POINT. dataset

N (per timestep)

timesteps

size

HACC (hiRes) HACC (lowRes) EXAALT FPM (hiRes) FPM (midRes) FPM (lowRes)

1,073,734,015 280,953,867 2,869,440 1,686,160 545,678 196,066

16 1 1 60 121 121

192.00 GB 3.14 GB 32.84 MB 1.13 GB 775.62 MB 271.50 MB

For each HACC (hiRes) timestep, the 2563 domain is partitioned across 64 MPI ranks using a hierarchical, interleaved decomposition. While the z and y dimensions follow standard linear partitioning (4 and 8 segments, respectively), the x-dimension employs a non-contiguous Fig. 2. Hierarchical spatial decom- strategy where each rank position of the 2563 HACC simulation volume across 64 MPI ranks. manages four disjoint subColor-coded blocks represent indi- intervals. As shown in Fig. 2, vidual rank assignments, with rank the domain is divided into 32 IDs annotated on visible faces. yz “super-blocks,” each split along x into two interleaved rank sets. This configuration creates spatially periodic “slabs” or “pencils” to facilitate long-range force calculations [46]. Hardware. Experiments were conducted on the National Energy Research Scientific Computing Center (NERSC) Perlmutter supercomputer [49]. GPU-based evaluations used NVIDIA A100 (40 GB HBM2) nodes with CUDA 12.4, while CPU benchmarks ran on 64-core AMD EPYC 7763 nodes. 63 31 59 63 27 55 31 59 23 51 63 27 55 19 47 23 51 15 43 31 59 11 39 63 27 55 19 47 23 51 15 43 59 7 11 39 55 19 47 51 15 43 7 31 35 11 39 47 27 7 43 23 35 3 39 3 19 30 34 15 35 26 2 11 3 22 7 34 35 2 18 29 3 33 14 25 10 2 34 33 1 21 6 34 17 1 32 28 2 13 33 24 0 9 1 20 32 5 33 16 0 1 32 12 8 0 4 32 0

B. Baseline comparison This section evaluates our proposed algorithm against state-of-the-art compressors, focusing on storage efficiency, throughput, and accuracy. Due to its substantial memory footprint exceeding a single GPU’s capacity, all HACC (hiRes) results were obtained using a distributed configuration of 16 nodes (4 GPUs per node) on a single timestep unless otherwise specified. Peak GPU memory consumption per rank remains around 1.144 GB across all tested configurations, well within the 40 GB A100 capacity. Experiments for all other datasets were conducted on a single GPU. In PGD, we apply the adaptive moment estimation (Adam) algorithm [44] to adaptively adjust learning rate for faster convergence. Adam is selected for its fast convergence and memory efficiency, which are critical for processing large-scale data on GPU. Following standard practice, we use Adam with α = 10−3 , β1 = 0.9, β2 = 0.999, ϵ = 10−8 ; results are robust to reasonable variation in α due to Adam’s adaptive gradient scaling [44]. The linking length b is a user-specified parameter; all experiments use η = 0.2 (cosmological standard).

Fig. 3. Compression ratio vs. maximum relative error of our method and baselines before and after correction on six datasets.

Fig. 4. Throughput of our correction method and baselines.

Storage Overhead. Fig. 3 illustrates the relationship between the compression ratio and the maximum relative error for all evaluated base compressors, both before and after our correction phase with full convergence loss tolerance Fig. 5. PSNR vs. BPP of our method under and baselines before and after correc- ϵL = 10−10 . Among the tion on EXAALT dataset. candidates, LCP yields the highest compression ratio across all error bounds, followed by Morton-order baseline and Draco. The storage overhead required to preserve FoF clusters increases as the relative error bound loosens, whereas the base compressor’s compression ratio improves with higher error tolerance. This competing dynamic creates an optimal point (or “peak”) in the effective compression ratio, as seen in the HACC (hiRes), HACC (lowRes), and EXAALT datasets. For example, the total compression ratio for LCP on the EXAALT data peaks at a relative error of approximately 6 × 10−4 . This suggests that if a baseline compression uses a coarser bound (e.g., 10−3 ), tightening the bound toward this peak value actually improves storage efficiency by reducing the subsequent correction overhead. Conversely, on FPM datasets, where vulnerable pairs are sparse, the storage overhead introduced by our correction is almost negligible, resulting in a monotonic increase in compression ratio without a discernible peak. The results on FPM datasets confirm that low violation density yields negligible overhead, and demonstrate applicability beyond particle physics to fluid dynamics. Rate distortion. We evaluate reconstruction accuracy versus storage cost by rate-distortion (RD) curves plotting PSNR against bits per particle (BPP) in Fig. 5. The reported BPP includes both the base compressed data and our edits, ensuring

a fair comparison. Our correction shifts RD curves upward or leaves them unchanged, with a substantial gap (more than 50 dB) for Draco. Since our method optimizes for FoF cluster membership rather than pointwise accuracy and may perturb particles away from their base-compressed positions, the observed RD improvement is a beneficial byproduct that demonstrates our correction reclaims accuracy lost during initial compression without compromising cluster fidelity. Throughput. The throughput of the 1500 { = 10-3 base compressors and our correction r ..)000 method with full convergence is illus� 500 trated in Fig. 4. Since HACC (hiRes) 0 dataset needs a distributed MPI con100 50 0 figuration while other datasets are pro2.0 L {= 10cessed on a single GPU, we provide 1.5 g, 1.0 a detailed scaling analysis for HACC � 0.5 (hiRes) in Section IV-D and focus here 0.0 on single-node performance. 100 50 0 num. iterations Our correction phase consistently Fig. 6. Tight loss Ltight outperforms the base compressors on vs. iteration count FPM datasets, achieving a speedup (limited to 100) for our of 6-25× over CPU-based baselines method correcting ZFPcompressed EXAALT and 15-30× over the GPU-accelerated data with relative cuSZp2. Conversely, for high-density ξ = 10−3 (top) and −4 datasets such as HACC (lowRes) and ξ = 10 (bottom). EXAALT, cuSZp2’s throughput exceeds our iterative refinement, as our method’s complexity scales with the density of vulnerable pairs. The correction runtime is primarily governed by the vulnerable pair count rather than the PGD iteration count alone; large vulnerable pair counts significantly increase per-iteration latency, whereas high iteration counts with small vulnerable pair counts remain computationally inexpensive. Table II illustrates this using ZFP-compressed EXAALT data. The execution times for ξ = 10−3 and ξ = 10−4 diverge significantly despite nearly identical iteration counts, due to a 22× difference in .c:

I

I

I

4

.c:

I

I

I

Fig. 7. Accuracy in FoF clustering results. (a) Matthews correlation coefficient (MCC) of link existence between particle pairs vs. bits per particle (BPP). (b) Halo mass functions (HMFs) and their pointwise relative errors for base compressors compared to our correction method at equivalent compression ratios. For x-values extending beyond the range of original data, relative error is calculated using linear interpolation of the original HMF.

vulnerable pair counts, confirming that throughput is sensitive to the number of violated pairs processed per iteration. Furthermore, we identify an approximately 0.25 s performance floor for spatial partitioning and vulnerable pair discovery, which persists even with zero violated pairs and iterations, serving as a verification step to ensure all particle pairs satisfy the prescribed linking-length constraints. TABLE II C ONVERGENCE OF THE CORRECTION ALGORITHM ON THE EXAALT DATASET. T IMES IN SECONDS ; ALL VIOLATIONS ARE ELIMINATED WITHIN REPORTED ITERATIONS . T HE QUANTIZATION STEP INTRODUCES NO NEW VIOLATIONS BECAUSE THE SAFETY MARGIN ϵq IN LTIGHT EXCEEDS THE MAXIMUM QUANTIZATION ERROR BY CONSTRUCTION . D OUBLE - PRECISION IS USED FOR DISTANCES TO ACCOMMODATE SMALL ERROR BOUNDS . rel. ξ

# vul. prs

# viol. prs

# iter

ttotal

tsetup

tP GD

10−3

41,445,603 1,876,427 185,048 18,508 1,808 200 19

305,855 37,959 4,783 313 34 0 0

105 102 4 2 1 0 0

2.01 0.870 0.295 0.280 0.262 0.257 0.255

0.253 0.261 0.246 0.250 0.253 0.257 0.255

1.757 0.609 0.049 0.030 0.009 0 0

10−4 10−5 10−6 10−7 10−8 10−9

Matthews correlation coefficient (MCC). Although our algorithm is proven to converge, the computational cost to reach full convergence may be prohibitive for datasets with high vulnerable-pair densities. For example, correcting ZFPcompressed HACC (lowRes) data with ξ = 10−3 requires 3,056 iterations to resolve 1.35 billion vulnerable pairs, taking 49.34 s, which is approximately twice as slow as ZFP’s compression process. However, as illustrated in Fig. 6, the loss decreases sharply in the initial steps before entering a long, slow convergence tail. Both factors, the high runtime cost of full convergence and the rapid early reduction in loss, suggest early termination as a viable heuristic. To evaluate the trade-off between execution time and clustering fidelity, we investigate the impact of a restricted iteration budget on FoF accuracy. We fix the maximum iteration count to 100 and plot MCC against bit-rate in BPP in Fig. 7 (a). Results for FPM datasets are omitted as the algorithm consistently converges within 100 iterations, yielding a trivial MCC of 1.0. Base compressors exhibit low MCC on the HACC (hiRes) dataset because they process ranks in isolation; this lack of inter-node coordination results in fractured halo structures at sub-domain boundaries. While our method may not reach strict convergence within the 100-iteration limit, it significantly outperforms the base compressors, demonstrating that even a truncated correction phase effectively restores the FoF clusters. Halo mass function (HMF). To validate the efficacy of the truncated correction process for downstream FoF analysis, we calculate the HMFs for results using a 100-iteration limit; as Table II shows, most configurations converge in far fewer iterations, and the cap is only approached under the most aggressive error bounds where residual loss is already negligible. Fig. 7 (b) compares the HMFs (B = 50 bins) derived from SZ3 and LCP, the base compressors with the highest compression ratios, against our correction results. The impact of our algorithm is most evident in the lower panel, which displays the pointwise relative error of the reconstructed HMFs. At the same compression ratios, our method consistently reduces the relative error across the whole HMFs, more accurately preserving the statistical distribution of the cosmological halos. C. Single-node performance Table III details per-kernel performance for our correction process on LCP-compressed EXAALT data (ξ = 10−4 , 100 iterations), involving 1.8M vulnerable pairs and 1.7M editable particles. The GPU achieves a 62× end-to-end speedup over the 64-core CPU baseline. All computational phases exhibit arithmetic intensity (AI) significantly below the roofline ridge points (A100: 9.75 F/B; EPYC: 24.5 F/B), confirming the algorithm is memory-bound in both cases. The three PGD kernels reach 43-49% of the A100’s peak HBM2 bandwidth (2 TB/s). The PGDUpdatePositions kernel assigns one GPU thread per editable particle, enabling embarrassingly parallel Adam updates and box-constrained projections. The PGDComputeLoss kernel uses sharedmemory binary-tree reduction to accumulate pair violations without global atomic traffic, sustaining near-peak bandwidth.

TABLE III S INGLE - NODE CUDA K ERNELS AND CPU F UNCTIONS P ERFORMANCE M ETRICS Kernel / Function

Calls

spatialPartitioning

1

findVulnerablePairs

2

PGDComputeLoss

10

PGDComputeGradients

100

PGDUpdatePositions

100

LosslesslyCompressEdits

1

total

1

Platform

Time

BW (GB/s)

BW Eff. (%)

AI (F/B)

GPU CPU GPU CPU GPU CPU GPU CPU GPU CPU GPU CPU GPU CPU

125 µs 1.75 s 256.9 ms 11.909 s 1.6 ms 1.664 s 15.7 ms 3.826 s 131.7 ms 3.220 s 4.3 ms 3.177 s 410 ms 25.5 s

354.5 1.2 46.1 1.2 971.7 30.4 984.2 22.2 872.6 22.4 43.1 1.2 351 8.9

17.7 0.6 2.3 0.6 48.6 14.8 49.2 10.8 43.6 10.9 2.2 0.6 17.6 4.3

2.25 0.005 6.92 0.038 0.41 0.032 0.17 0.019 0.12 0.127 0.05 0.003 0.69 0.057

The full GPU pipeline’s effective bandwidth is 351 GB/s (17.6% of peak), with the reduction from kernel-level peaks due to the structural irregularity of the spatial search. The findVulnerablePairs kernel is the primary bottleneck (63% of GPU time, 47% of CPU time) and reflects the complexity of distance-constraint enumeration in irregular particle data. With an AI of 6.92 F/B but only 2.3% bandwidth utilization, the kernel is latency-bound due to address-dependent hash-table probe. Our GPU implementation mitigates this through three design choices: (1) CUB-based spatial sorting and adaptive grid coarsening to ensure spatially contiguous thread access; (2) a two-pass atomic vulnerable pair counting strategy for exact memory allocation; and (3) local thread-level buffering to minimize global atomic contention. Furthermore, a custom flat hash table using atomicCAS replaces the CPU’s std::unordered_map, providing O(1) lookups. These co-design optimizations yield a 46× speedup and increase bandwidth efficiency from 10.8% (CPU) to 49.2% (GPU).

(a)

(b)

Fig. 8. Strong scaling analysis for the GPU implementation on the HACC (hiRes) dataset. (a) Wall-clock time (blue line) and parallel efficiency (texts) relative to the 16-process baseline; the dashed line represents ideal linear scaling. (b) Analysis of system imbalance across three metrics: execution time per rank, local particle distribution, and the ghost-to-local particle ratio.

D. Inter-node distributed scaling This section evaluates CPU-MPI and GPU-MPI performance to illustrate how distributed-memory coordination scales for exascale-class workloads. Strong scaling. Fig. 8 and Table IV present a strong scaling analysis of our correction algorithm on SZ3-compressed HACC (hiRes) dataset with relative ξ = 10−6 in overall and breakdown senses, respectively. To evaluate 16-256 processes

speedup 14,000× 46.4× 1,040× 244× 24.4× 739× 62×

while maintaining HACC’s geometric logic as in in-situ compression and correction, we aggregate data by merging ranks with adjacent indices into larger units (e.g., grouping ranks 0–3 into a single process for a 16-process run). Conversely, for 128- and 256-process configurations requiring higher granularity, the data within each rank is subdivided along the xaxis. Since HACC’s spatial decomposition changes with rank count (adjacent ranks merged for fewer nodes, sub-divided for more), scaling results reflect both parallel efficiency and decomposition effects; a fixed-decomposition strong scaling test would require re-partitioning the simulation data, which is outside the scope of this work. Correctness is preserved across √ all scales: the ghost zone width δ = b + 2 3ξ guarantees every boundary-crossing halo and its PGD gradient are fully replicated, making the results invariant to process count. GPU acceleration and overall speedup. Table IV shows that the GPU-implementation achieves end-to-end speedups of 62.8×, 56.6×, and 55.2× over the CPU at 16, 64, and 256 processes, respectively. The modest decline in speedup ratio from 62.8× to 55.2× as scale increases reflects Amdahl’s law [50]: at 16 ranks the GPU eliminates nearly all local compute cost, but at 256 ranks the communication overhead constitutes a larger fraction of total runtime. Strong scaling efficiency. Both implementations scale well across the 16× process increase. The CPU achieves 82.0% from 16 to 256 ranks. As shown in Fig. 8 (a), GPU efficiency degrades faster than that of the CPU because the fast GPU compute makes the slowly-scaling communication and I/O the primary bottlenecks of the total runtime. From 16 and 256 ranks, local compute remains stable (∼30% of runtime), but communication cost grows from 41.2% to 44.0% and I/O time grows from 14.9% to 22.4%. This shift toward a communication- and I/O-bound regime explains why the CPU, dominated by local compute, exhibits superior scaling. Convergence check calls MPI_Allreduce to reduce a single float per iteration; the reported time (0.09–0.14 s over 100 iterations) is almost entirely barrier synchronization. Imbalance exceeds 100% at every scale on both GPU and CPU, meaning the slowest rank takes over 2× the mean, which is driven by uneven local and ghost particle distribution

TABLE IV P ERFORMANCE BREAKDOWN AND LOAD IMBALANCE ACROSS SCALES FOR CPU AND GPU MPI IMPLEMENTATIONS ON HACC ( HI R ES ). T IMES ARE MAXIMUM SECONDS ACROSS RANKS ; IMBALANCE IS (tmax − tavg )/tavg . CPU Time (sec)

CPU Imbalance

GPU Time (sec)

GPU Imbalance

Group

Phase

16p

64p

256p

16p

64p

256p

16p

64p

256p

16p

64p

256p

I/O

Read original Read decompressed Write edits

0.79 0.74 0.02

0.22 0.21 0.01

0.17 0.05 0.00

7.6% 7.4% 17.9%

16.8% 13.6% 35.5%

217.8% 41.3% 81.3%

0.39 0.74 0.02

0.12 0.23 0.01

0.07 0.07 0.01

14.7% 6.27% 6.8%

32.5% 311.3% 17.0%

528.6% 50.8% 54.6%

Comm.

BBox exch. Ghost exch. Converg. check

0.36 4.89 0.14

0.40 1.63 0.14

0.46 0.68 0.08

0.02% 10.4% 96.1%

0.3% 6.2% 135.6%

0.9% 13.1% 227.1%

0.06 3.04 0.09

0.08 0.334 0.10

0.12 0.115 0.06

95.1% 17.5% 116.2%

120.2% 16.6% 116.6%

152.5% 106.7% 117.9%

Local Compute

Grid partition VP detection PGD optim. Huffman + ZSTD

58.22 251.9 78.66 89.98

13.27 59.28 28.52 20.23

3.15 16.15 11.10 5.16

2.3% 9.7% 53.8% 7.5%

10.3% 20.9% 52.8% 17.8%

30.8% 53.9% 37.3% 48.5%

0.126 0.732 0.85 0.567

0.01 0.18 0.636 0.217

0.005 0.061 0.12 0.02

15.0% 24.8% 7.4% 5.3%

34.3% 31.9% 21.8% 19.2%

96.5% 39.8% 26.9% 58.3%

H2D / D2H

—

—

—

—

—

—

1.125

0.273

0.019

8.0%

20.6%

129.6%

Total Wall Clock

485.7

124.0

37.0

—

—

—

—

—

—

End-to-end GPU Speedup

—

(Fig. 8 (b)). Overlapping the reduction with local computation via MPI_Iallreduce could recover this idle time. I/O imbalance. At high ranks, Lustre filesystem contention [51] becomes the primary scaling bottleneck. GPU read-input imbalance reaches 528.6% at 256 ranks, with the slowest rank taking 6.3× the mean, stalling the entire system at subsequent barriers and limiting efficiency beyond this scale.

(a)

(b)

Fig. 9. Weak scaling execution time breakdown and throughput for the GPU implementation. (a) Wall-clock time breakdown across scales. Each stacked bar shows the maximum time over all ranks. Percentages on the bars are overall weak scaling efficiency. (b) The aggregate throughput. The dashed line represents ideal linear scaling. TABLE V P ERFORMANCE AND W ORKLOAD M ETRICS FOR W EAK S CALING . Proc.

Avg. Particles per Rank

Data Imbalance

Time Imbalance

64 128 256 512 1,024

16,777,117 16,777,115 16,777,109 16,777,101 16,777,083

12.22% 12.93% 14.54% 16.56% 24.69%

1.154% 1.337% 1.406% 1.638% 1.568%

Weak scaling. To simulate exascale-class workloads, we concatenate HACC (hiRes) timesteps along the x-axis, maintaining a constant per-rank workload. The problem scales from 16 nodes (64 ranks, 1 dataset) up to 256 nodes (1,024 ranks, 16 datasets), increasing the simulation domain from 2563 to 4, 096 × 256 × 256. While not physically equivalent to larger simulation volumes, this data construction of increasing size

7.74

2.19

0.67

62.8×

56.6×

55.2×

—

isolates how correction cost scales with N , independent of variations in clustering structure or simulation parameters. Fig. 9 shows the execution time breakdown, weak scaling efficiency, and aggregate throughput as the problem scales from 64 to 1,024 ranks. Our method maintains a efficiency of 72.6% at the 1,024-rank scale (Fig. 9 (a)). While local compute remains stable, wall-clock growth is driven by communication overhead, which increases with rank count. This is characteristic of 3D decompositions in cosmological simulations, where high concurrency elevates the surface-to-volume ratio and communication pressure. Despite this, the throughput scales near-ideally, reaching ∼62 GB/s at peak concurrency. Table V summarizes the per-rank workload distribution across the scaling range. The average particle count per rank remains consistent at ∼ 1.68 × 107 , verifying the integrity of the weak-scaling problem construction. Although the data imbalance grows from 12.22% (64 ranks) to 24.69% (1,024 ranks), the time imbalance remains stably below 1.7% throughout, peaking at only 1.638% at 512-rank and even slightly regressing at the 1,024-rank scale. This performance decoupling suggests that our GPU kernels effectively mask data-level skew through high occupancy and efficient memory access, preventing workload imbalance from stalling the critical path. V. C ONCLUSION We presented a high-performance algorithm that preserves FoF cluster membership under lossy compression by formulating connectivity correction as a constrained optimization problem solved via our projected gradient descent. Our implementation leverages GPU parallelism and MPI coordination to handle exascale-class workloads, and experiments on NERSC’s Perlmutter demonstrate effective cluster membership restoration across cosmological, molecular dynamics, and fluid dynamics datasets with negligible impact on compression ratio. While convergence is guaranteed when violations consist solely of broken links, the common case at moderate error bounds, aggressive bounds increase costs as complex configurations may require hundreds of iterations.

Limitations. The correction method is specific to the linking length b specified at compression time; re-analysis with a different b requires reapplying the correction to the base decompressed data. At aggressive error bounds (ξ ≥ 10−3 ), correction cost can exceed base compression time when vulnerable pair density is high. Additionally, the current implementation assumes a non-periodic domain and does not support 6D phase-space finders like Rockstar. Future works. First, we will extend this correction algorithm to time-evolving datasets, leveraging temporal coherence to reduce the edit footprint. Second, we aim to generalize the loss function to phase-space halo finders such as Rockstar. Third, we plan to explore integrating the correction directly into the base compressor’s transform stage to eliminate the secondary correction log. ACKNOWLEDGMENT The material was supported by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research (ASCR), under contracts DE-AC02-06CH11357 and DE-SC0025677. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a Department of Energy User Facility using NERSC award DDR-ERCAP 0034457. R EFERENCES [1] S. Habib, A. Pope, H. Finkel, N. Frontiere, K. Heitmann, D. Daniel, P. Fasel, V. Morozov, G. Zagaris, T. Peterka et al., “HACC: Simulating sky surveys on state-of-the-art supercomputing architectures,” New Astronomy, vol. 42, pp. 49–65, 2016. [Online]. Available: https: //doi.org/10.1016/j.newast.2015.06.003 [2] N. Frontiere, J. D. Emberson, M. Buehlmann, E. M. Rangel, S. Habib, K. Heitmann, P. Larsen, V. Morozov, A. Pope, C.-A. Faucher-Giguère et al., “Cosmological hydrodynamics at exascale: A trillion-particle leap in capability,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2025, pp. 25–35. [Online]. Available: https://doi.org/10.1145/3712285.3771786 [3] N. Tchipev, S. Seckler, M. Heinen, J. Vrabec, F. Gratl, M. Horsch, M. Bernreuther, C. W. Glass, C. Niethammer, N. Hammer et al., “TweTriS: Twenty trillion-atom simulation,” The International Journal of High Performance Computing Applications, vol. 33, no. 5, pp. 838–854, 2019. [Online]. Available: https://doi.org/10.1177/1094342018819741 [4] T. Xie, A. France-Lanord, Y. Wang, J. Lopez, M. A. Stolberg, M. Hill, G. M. Leverick, R. Gomez-Bombarelli, J. A. Johnson, Y. Shao-Horn et al., “Accelerating amorphous polymer electrolyte screening by learning to reduce errors in molecular dynamics simulated properties,” Nature communications, vol. 13, no. 1, p. 3415, 2022. [Online]. Available: https://doi.org/10.1038/s41467-022-30994-1 [5] Z. Jian, Y. Chen, S. Xiao, L. Wang, X. Li, K. Wang, H. Deng, and W. Hu, “Shock-induced plasticity and phase transformation in single crystal magnesium: an interatomic potential and non-equilibrium molecular dynamics simulations,” Journal of Physics: Condensed Matter, vol. 34, no. 11, p. 115401, 2022. [Online]. Available: https://doi.org/10.1088/1361-648X/ac443e [6] K. Zhao, S. Di, M. Dmitriev, T.-L. D. Tonellot, Z. Chen, and F. Cappello, “Optimizing error-bounded lossy compression for scientific data by dynamic spline interpolation,” in Proceedings of 2021 IEEE 37th International Conference on Data Engineering (ICDE). IEEE, 2021, pp. 1643–1654. [Online]. Available: https: //doi.org/10.1109/ICDE51399.2021.00145

[7] S. Jin, D. Tao, H. Tang, S. Di, S. Byna, Z. Lukic, and F. Cappello, “Accelerating parallel write via deeply integrating predictive lossy compression with HDF5,” in SC22: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2022, pp. 1–15. [Online]. Available: https: //doi.org/10.1109/SC41404.2022.00066 [8] C. Ren, R. Underwood, S. Di, E. Kutay, Z. Lukic, A. Yener, F. Cappello, and H. Guo, “FFCz: Fast Fourier Correction for SpectrumPreserving Lossy Compression of Scientific Data,” arXiv preprint arXiv:2601.01596, 2026. [Online]. Available: https://doi.org/10.48550/ arXiv.2601.01596 [9] L. Zhang, R. Li, C. Ren, S. Di, J. Liu, J. Huang, R. Underwood, P. Grosset, D. Tao, X. Liang et al., “LCP: Enhancing Scientific Data Management with Lossy Compression for Particles,” Proceedings of the ACM on Management of Data, vol. 3, no. 1, pp. 1–27, 2025. [Online]. Available: https://doi.org/10.1145/3709700 [10] X. Liang, H. Guo, S. Di, F. Cappello, M. Raj, C. Liu, K. Ono, Z. Chen, and T. Peterka, “Toward feature-preserving 2D and 3D vector field compression,” in PacificVis, 2020, pp. 81–90. [Online]. Available: https://doi.org/10.1109/PacificVis48177.2020.6431 [11] X. Liang, S. Di, F. Cappello, M. Raj, C. Liu, K. Ono, Z. Chen, T. Peterka, and H. Guo, “Toward feature-preserving vector field compression,” IEEE Transactions on Visualization and Computer Graphics, vol. 29, no. 12, pp. 5434–5450, 2022. [Online]. Available: https://doi.org/10.1109/TVCG.2022.3214821 [12] X. Liang, K. Zhao, S. Di, S. Li, R. Underwood, A. M. Gok, J. Tian, J. Deng, J. C. Calhoun, D. Tao et al., “SZ3: A modular framework for composing prediction-based error-bounded lossy compressors,” IEEE Transactions on Big Data, vol. 9, no. 2, pp. 485–498, 2022. [Online]. Available: https://doi.org/10.1109/TBDATA.2022.3201176 [13] P. Lindstrom, “Fixed-rate compressed floating-point arrays,” IEEE transactions on visualization and computer graphics, vol. 20, no. 12, pp. 2674–2683, 2014. [Online]. Available: https://doi.org/10.1109/ TVCG.2014.2346458 [14] Google, “Google Draco,” https://github.com/google/draco, 2024, accessed Feb. 08, 2026. [15] S. More, A. V. Kravtsov, N. Dalal, and S. Gottlöber, “The overdensity and masses of the friends-of-friends halos and universality of halo mass function,” The Astrophysical Journal Supplement Series, vol. 195, no. 1, p. 4, 2011. [Online]. Available: https: //doi.org/10.1088/0067-0049/195/1/4 [16] F. Rodriguez and M. Merchán, “Combining friend-of-friend and halo-based algorithms for the identification of galaxy groups,” Astronomy & Astrophysics, vol. 636, p. A61, 2020. [Online]. Available: https://doi.org/10.1051/0004-6361/201937423 [17] R. Monchaux, M. Bourgoin, and A. Cartellier, “Preferential concentration of heavy particles: a Voronoı̈ analysis,” Physics of Fluids, vol. 22, no. 10, 2010. [Online]. Available: https://doi.org/10.1063/1.3489987 [18] J. R. West, T. Maurel-Oujia, K. Matsuda, K. Schneider, S. S. Jain, and K. Maeda, “Clustering, rotation, and swirl of inertial particles in turbulent channel flow,” International Journal of Multiphase Flow, vol. 174, p. 104764, 2024. [Online]. Available: https: //doi.org/10.1016/j.ijmultiphaseflow.2024.104764 [19] A. Colanera, J. M. Reumschüssel, J. P. Beuth, M. Chiatto, L. De Luca, and K. Oberleithner, “Extended cluster-based network modeling for coherent structures in turbulent flows,” Theoretical and Computational Fluid Dynamics, vol. 39, no. 1, p. 1, 2025. [Online]. Available: https://doi.org/10.1007/s00162-024-00723-z [20] A. France-Lanord and J. C. Grossman, “Correlations from ion pairing and the Nernst-Einstein equation,” Physical review letters, vol. 122, no. 13, p. 136001, 2019. [Online]. Available: https: //doi.org/10.1103/PhysRevLett.122.136001 [21] C. Li, X. Fu, W. Zhong, and J. Liu, “Dissipative particle dynamics simulations of a protein-directed self-assembly of nanoparticles,” ACS omega, vol. 4, no. 6, pp. 10 216–10 224, 2019. [Online]. Available: https://doi.org/10.1021/acsomega.9b01078 [22] F. Cappello, M. Ainsworth, J. Bessac, M. Burtscher, J. Y. Choi, E. Constantinescu, S. Di, H. Guo, P. Lindstrom, and O. Tugluk, “SDRBench: Scientific Data Reduction Benchmarks,” https://sdrbench. github.io/, 2020, online; accessed Jul. 08, 2025. [23] Z. Lukić, K. Heitmann, S. Habib, S. Bashinsky, and P. M. Ricker, “The halo mass function: High-redshift evolution and universality,” The

Astrophysical Journal, vol. 671, no. 2, pp. 1160–1181, 2007. [Online]. Available: https://doi.org/10.1086/523083 [24] M. S. Warren, K. Abazajian, D. E. Holz, and L. Teodoro, “Precision determination of the mass function of dark matter halos,” The Astrophysical Journal, vol. 646, no. 2, pp. 881–885, 2006. [Online]. Available: https://doi.org/10.1086/504962 [25] V. Springel, S. D. White, A. Jenkins, C. S. Frenk, N. Yoshida, L. Gao, J. Navarro, R. Thacker, D. Croton, J. Helly et al., “Simulations of the formation, evolution and clustering of galaxies and quasars,” nature, vol. 435, no. 7042, pp. 629–636, 2005. [Online]. Available: https://doi.org/10.1038/nature03597 [26] S. Di, J. Liu, K. Zhao, X. Liang, R. Underwood, Z. Zhang, M. Shah, Y. Huang, J. Huang, X. Yu et al., “A survey on errorbounded lossy compression for scientific datasets,” ACM computing surveys, vol. 57, no. 11, pp. 1–38, 2025. [Online]. Available: https://doi.org/10.1145/3733104 [27] Y. Huang, S. Di, G. Li, and F. Cappello, “cuSZp2: A GPU lossy compressor with extreme throughput and optimized compression ratio,” in SC24: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2024, pp. 1–18. [Online]. Available: https://doi.org/10.1109/SC41406.2024.00021 [28] L. Ibarria, P. Lindstrom, J. Rossignac, and A. Szymczak, “Outof-core compression and decompression of large N-dimensional scalar fields,” in Computer Graphics Forum, vol. 22, no. 3. Wiley Online Library, 2003, pp. 343–348. [Online]. Available: https://doi.org/10.1111/1467-8659.00681 [29] MPEG Group, “MPEG G-PCC,” https://github.com/MPEGGroup/ mpeg-pcc-tmc13, 2024. [30] ——, “MPEG V-PCC,” https://github.com/MPEGGroup/ mpeg-pcc-tmc2, 2024. [31] K. Zhao, S. Di, D. Perez, X. Liang, Z. Chen, and F. Cappello, “MDZ: An efficient error-bounded lossy compressor for molecular dynamics,” in 2022 IEEE 38th International Conference on Data Engineering (ICDE). IEEE, 2022, pp. 27–40. [Online]. Available: https://doi.org/10.1109/ICDE53745.2022.00007 [32] J. Liu, P. Jiao, K. Zhao, X. Liang, S. Di, and F. Cappello, “QPET: A versatile and portable quantity-of-interest-preservation framework for error-bounded lossy compression,” Proceedings of the VLDB Endowment, vol. 18, no. 8, pp. 2440–2453, 2025. [Online]. Available: https://doi.org/10.14778/3742728.3742739 [33] Y. Li, X. Liang, B. Wang, Y. Qiu, L. Yan, and H. Guo, “MSz: An efficient parallel algorithm for correcting morse-smale segmentations in error-bounded lossy compressors,” IEEE Transactions on Visualization and Computer Graphics, vol. 31, no. 1, pp. 130–140, 2024. [Online]. Available: https://doi.org/10.1109/TVCG.2024.3456337 [34] N. Gorski, X. Liang, H. Guo, L. Yan, and B. Wang, “A general framework for augmenting lossy compressors with topological guarantees,” IEEE Transactions on Visualization and Computer Graphics, 2025. [Online]. Available: https://doi.org/10.1109/TVCG. 2025.3567054 [35] F. Roy, V. R. Bouillot, and Y. Rasera, “pFoF: a highly scalable halo-finder for large cosmological data sets,” Astronomy & Astrophysics, vol. 564, p. A13, 2014. [Online]. Available: https://doi.org/10.1051/ 0004-6361/201322555 [36] K. Heitmann, T. D. Uram, H. Finkel, N. Frontiere, S. Habib, A. Pope, E. Rangel, J. Hollowed, D. Korytov, P. Larsen et al., “HACC cosmological simulations: First data release,” The Astrophysical Journal Supplement Series, vol. 244, no. 1, p. 17, 2019. [Online]. Available: https://doi.org/10.3847/1538-4365/ab3724 [37] V. Springel, N. Yoshida, and S. D. White, “GADGET: a code for collisionless and gasdynamical cosmological simulations,” New Astronomy, vol. 6, no. 2, pp. 79–117, 2001. [Online]. Available: https://doi.org/10.1016/S1384-1076(01)00042-2 [38] D. S. Reed, R. Bower, C. S. Frenk, A. Jenkins, and T. Theuns, “The halo mass function from the dark ages through the present day,” Monthly Notices of the Royal Astronomical Society, vol. 374, no. 1, pp. 2–15, 2007. [Online]. Available: https://doi.org/10.1111/j.1365-2966. 2006.11204.x [39] J. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. Warren, G. Yepes, S. Gottlöber, and D. E. Holz, “Toward a halo mass function for precision cosmology: The limits of universality,” The Astrophysical Journal, vol. 688, no. 2, pp. 709–728, 2008. [Online]. Available: https://doi.org/10.1086/591439

[40] A. R. Zentner, “The excursion set theory of halo mass functions, halo clustering, and halo growth,” International Journal of Modern Physics D, vol. 16, no. 05, pp. 763–815, 2007. [Online]. Available: https://doi.org/10.1142/S0218271807010511 [41] D. A. Huffman, “A method for the construction of minimum-redundancy codes,” Proceedings of the IRE, vol. 40, no. 9, pp. 1098–1101, 1952. [Online]. Available: https://doi.org/10.1109/JRPROC.1952.273898 [42] Y. Collet, “Zstandard (ZSTD),” https://github.com/facebook/zstd, 2015, version 1.5.7, Online; accessed Jun. 15, 2025. [43] A. Beck, First-order methods in optimization. SIAM, 2017. [Online]. Available: https://doi.org/10.1137/1.9781611974997 [44] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014. [Online]. Available: https: //doi.org/10.48550/arXiv.1412.6980 [45] NVIDIA Corporation, “CUDA Unbound (CUB),” https://docs.nvidia. com/cuda/cub/index.html, 2024, accessed Jul. 08, 2025. [46] S. Habib, V. Morozov, N. Frontiere, H. Finkel, A. Pope, and K. Heitmann, “HACC: Extreme scaling and performance across diverse architectures,” in Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, 2013, pp. 1–10. [Online]. Available: https://doi.org/10.1145/2503210.2504566 [47] Exascale Computing Project, “EXAALT: Exascale Atomistics for Accuracy, Length, and Time,” https://www.exascaleproject.org/ research-project/exaalt/, 2021, accessed: 2024-05-20. [48] S. C. . Data, “Finite Pointset Method (FPM) Viscous Fingers dataset,” https://cloud.sdsc.edu/v1/AUTH sciviscontest/2016/README. html, 2016, accessed: 2026-03-27. [49] N. E. R. S. C. C. (NERSC), “Perlmutter architecture,” https://docs.nersc. gov/systems/perlmutter/architecture/, 2015, accessed: 2026-03-20. [50] J. L. Gustafson, “Reevaluating amdahl’s law,” Communications of the ACM, vol. 31, no. 5, pp. 532–533, 1988. [Online]. Available: https://doi.org/10.1145/42411.42415 [51] W. Yu, J. Vetter, R. S. Canon, and S. Jiang, “Exploiting Lustre file joining for effective collective IO,” in Seventh IEEE International Symposium on Cluster Computing and the Grid (CCGrid’07). IEEE, 2007, pp. 267–274. [Online]. Available: https://doi.org/10.1109/ CCGRID.2007.51