arXiv:2609.04573v1 [cs.DC] 4 Sep 2026
Improving Progressive Compression with Adaptive Interpolation and Coefficient Decomposition Wenbo Li
Xuan Wu
Qian Gong
Oregon State University Corvallis, USA [email protected]
Oregon State University Corvallis, USA [email protected]
Oak Ridge National Laboratory Oak Ridge, USA [email protected]
Pu Jiao
Jieyang Chen
Qing Liu
University of Kentucky Lexington, USA [email protected]
University of Oregon Eugene, USA [email protected]
New Jersey Institute of Technology Newark, USA [email protected]
Norbert Podhorszki
Scott Klasky
Xin Liang
Oak Ridge National Laboratory Oak Ridge, USA [email protected]
Oak Ridge National Laboratory Oak Ridge, USA [email protected]
Oregon State University Corvallis, USA [email protected]
Abstract—Exascale simulations generate data far faster than it can be stored or analyzed, making efficient data reduction essential. Error-controlled lossy compression offers high compression ratios under user-specified error bounds, but the target tolerance must be fixed at compression time. Progressive compression relaxes this restriction, yet existing methods still rely on fixed refactoring strategies and do not fully exploit correlations among decomposed coefficients, limiting the efficiency of progressive retrieval. In this work, we present an adaptive progressive compression framework that improves retrieval efficiency for two common targets, namely error-bound and peak Signal-toNoise ratios. Our contributions are fourfold. (1) We propose to leverage two complementary interpolation schemes for adaptive progressive compression toward different targets, and we optimize them to achieve high efficiency. (2) We propose coefficient decomposition, a novel method that exploits the commonly overlooked spatial correlations among decorrelated data, which further improves the efficiency. (3) We develop the adaptive progressive compression workflow with automatic selection of the best-fit refactoring pipeline and tailored optimizations. (4) We evaluate the proposed framework on five real-world scientific datasets against three state-of-the-art progressive compressors. Experimental results demonstrate that the proposed framework improves the compression ratio by up to 42.3% under the same requested error tolerance and up to 92.5% at the same PSNR, compared with the best-performing existing methods. When transferring 512 GB of scientific data to remote sites, the framework delivers up to 1.26× speedup in the end-to-end data transfer performance. Furthermore, our method achieves the highest visualization quality while retrieving the least amount of data from storage.
This manuscript has been authored in part by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a non-exclusive, paid up, irrevocable, worldwide license to publish or reproduce the published form of the manuscript, or allow others to do so, for U.S. Government purposes. The DOE will provide public access to these results in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).
SC26, November 15-20, 2026, Chicago, Illinois, USA 979-8-3195-4789-7/26/$31.00 ©2026 IEEE
Index Terms—High performance computing, data compression, data processing, scientific computing
I. I NTRODUCTION Modern exascale simulations can generate data much faster than it can be stored, transferred, or analyzed. For example, a recent direct numerical simulation of isotropic turbulence on Frontier [1] produced about 0.5 PB per snapshot and would generate up to ∼38 TB/s if every time step were retained, far exceeding the 1–2 TB/s bandwidth of modern parallel file systems [2]. This growing gap makes efficient data reduction essential for computer and computational science. Error-controlled lossy compression [3]–[7] is a widely used approach to mitigate this bottleneck. Compressors such as SZ [4], [8]–[10], ZFP [6], and MGARD [7], [11]–[13] provide much higher compression ratios than lossless methods while guaranteeing that reconstruction error stays within a user-specified tolerance, and they have been integrated into scientific data management libraries such as HDF5 [14] and ADIOS-2 [15]. However, a key limitation of these singleerror-bounded lossy compressors is that users must specify an error bound before storing the data. Since the actual error tolerance required by various downstream analyses is often unknown a priori, and details lost during compression cannot be recovered, users tend to set conservative error bounds, thereby sacrificing compression ratio to preserve more information and reducing the overall effectiveness of lossy compression. Progressive data compression [16]–[22] tackles the above limitation by providing a more flexible reconstruction strategy. The idea originates from progressive image coding standards such as JPEG [19] and JPEG2000 [20], where data are refactored once into a near-lossless representation and then incrementally reconstructed until the desired accuracy is reached
during retrieval. These standards have long been used in modern browsers for progressive image rendering. However, they are not suitable for scientific applications because they cannot guarantee numerical error bounds. To accommodate scientific use cases, PMGARD [16] combines the MGARD decomposition theory [12] with bitplane encoding to offer guaranteed error control during progressive retrieval. Nonetheless, its efficiency is limited by the multilinear decomposition it adopts and by its greedy retrieval strategy, both of which leave room for improvement. Recently, IPComp [17] achieves better efficiency by employing cubic spline interpolation and a dynamic programming-based retrieval strategy. However, both PMGARD and IPComp are built around fixed interpolation pipelines, which limit their ability to adapt to different target error metrics and data characteristics. In addition, existing progressive compressors typically encode decomposed coefficients independently, overlooking their correlations and thereby leading to suboptimal efficiency. In this work, we design and develop an adaptive refactoring framework that improves retrieval efficiency for the two most commonly used types of target error tolerance. We further propose a novel method, namely coefficient decomposition, to improve efficiency by exploiting the spatial correlations within decomposed coefficients. We also develop an adaptive progressive compression workflow with online parameter tuning and tailored performance optimizations. Our contributions are summarized as follows: We leverage and optimize two complementary interpolation schemes for adaptive data decomposition in progressive compression, which offers high flexibility to achieve high efficiency towards diverse targets. • We propose a novel coefficient decomposition method that exploits spatial correlations within the decorrelated data after interpolation. This further improves retrieval efficiency with negligible metadata overhead. • We develop a fully adaptive progressive compression workflow that uses a sampling-based method to identify the best-fit refactoring pipeline on the fly and several tailored optimizations to mitigate overhead. • We rigorously evaluate our method against state-of-the-art progressive compression methods through comprehensive experiments on five real-world datasets. Experimental results demonstrate that our method achieves competitive performance, the best retrieval efficiency, and good scalability, showing up to 42.3% and 92.5% improvement in compression ratio under the same error tolerance and PSNR, respectively, 1.26× speedup in end-to-end transfer of 512 GB data, and the highest visualization quality with the least retrieved data size. •
The rest of the paper is organized as follows. Section II discusses the related works. Section III formulates the research problem and provides an overview of the proposed framework. Section IV introduces the adopted interpolation schemes. Section V describes the proposed coefficient decomposition in detail. Section VI summarizes the implementation and
optimization. Section VII presents the experimental evaluation with state-of-the-art methods and real-world datasets. Section VIII concludes the research with a vision for future work. II. R ELATED W ORK In this section, we recap related works on scientific lossy compressors and progressive compressors. A. Error-controlled Lossy Compressors Error-controlled lossy compressors target floating-point scientific data, where bounded numerical fidelity matters more than exact bitwise reconstruction. Unlike general-purpose lossless compressors such as GZIP [23], ZSTD [24], and BLOSC [25], they model spatial smoothness, multiscale structure, and value correlation in the data. Existing error-controlled compressors can be broadly grouped by how they decorrelate data before entropy coding. Prediction-based methods estimate each value from neighboring samples and then compress the residual. The SZ family [4], [8]–[10] is the best-known representative in this category. Early versions of SZ center on a Lorenzo predictor [26], while later versions incorporate richer models such as regression and interpolation-based predictors [8], [9]. After prediction, the residuals are quantized with a linear-scaling quantizer [4] and compressed through entropy encoding, such as Huffman coding [27], together with lossless coders like ZSTD [24]. Transform-based compressors instead map the data into another basis in which a small set of coefficients captures most of the data content. ZFP [6] is a representative design that partitions data into small blocks, converts values to fixed-point form under a common exponent, applies a near-orthogonal transform, and then uses embedded coding to emit only as many transformed bits as needed to satisfy requested accuracy. MGARD [7], [11], [12] provides another important line of work. It combines multilevel decomposition with finite-element analysis and derives error control from a mathematically grounded hierarchy. Despite their methodological differences, the prevailing usage model of these compressors remains the same, where compression is performed against a single user-specified tolerance, and the compressed representation is tailored to satisfy that tolerance alone. In practice, however, the same dataset often serves multiple downstream analyses with varying precision requirements. Since lost information cannot be recovered after decompression, users are compelled to select a conservatively tight error bound that satisfies the most demanding analysis, leading to unnecessarily low compression ratios and diminished benefits in data storage and transfer. B. Data Refactoring and Progressive Retrieval Progressive compression studies not only how data are compacted, but also how compressed information is organized so that reconstruction can proceed progressively. In this setting, a compressed representation is usually divided into ordered components such as bitplanes and hierarchical levels, so that it can be used to incrementally reconstruct the data to any
desired accuracy during retrieval. This write-once, retrieveprogressively strategy is particularly well-suited to scientific data management, where a dataset is generated once and subsequently consumed by numerous analyses with different fidelity requirements. The roots of this idea can be traced to progressive image coding methods such as JPEG and JPEG2000 [19], [20], which organize encoded information so that visual quality improves gradually as more bits arrive. Their influence on scientific compression lies mainly in progressive transmission and embedded coding concepts. PMGARD [16] is among the first compressors to bring this progressive retrieval paradigm into scientific error-bounded compression. Built on the MGARD theory [11], [12], it organizes multilevel coefficients into bitplane streams that can be decoded in stages while preserving strict error guarantees for both L2 and L∞ targets. However, PMGARD relies on linear interpolation during decomposition, which limits its ability to capture data correlations compared to higher-order interpolation schemes. Its greedy retrieval strategy also does not minimize retrieval volume for a given tolerance. Additionally, although PMGARD supports different target error metrics through per-bit and negabinary encodings, these changes are largely confined to the encoding phase, while the underlying decorrelation method remains the same. Magri et al. [28] explore a simpler construction that builds a progressive representation by repeatedly invoking an existing error-controlled lossy compressor with a sequence of decreasing tolerances. At a high level, the original data is first compressed at a coarse tolerance, then finer residuals are compressed successively to form additional refinement levels. While straightforward to implement, its drawback is that retrieval becomes accumulation-heavy, since reconstructing a fine level requires decoding all preceding levels and summing their contributions, causing the retrieval cost to grow with the number of refinement stages. IPComp [17] advances the state of the art by combining the interpolation-based decomposition from SZ3 [9] with a dynamic-programming retrieval strategy. Compared with PMGARD, its cubic spline interpolation captures stronger local correlation, and its retrieval algorithm is better aligned with minimizing retrieved data volume. Nevertheless, notable limitations remain. First, it quantizes the data with a preset error bound before decorrelation, thereby reducing the achievable precision and degrading decorrelation accuracy. Second, it leverages a fixed pipeline that cannot adapt to different error metrics. Last but not least, it overlooks the correlation in the decorrelated coefficients as existing work, which limits the overall efficiency. These shortcomings motivate the present work, in which we design an adaptive refactoring framework that addresses the above limitations through adaptive interpolation and encoding, a novel coefficient decomposition method that exploits spatial correlations within decomposed data, and a systematic design that delivers an adaptive progressive compression pipeline with tailored performance optimizations.
III. OVERVIEW In this section, we formulate the research problem and provide an overview of the proposed framework. A. Problem Formulation We define the error-controlled progressive method, and then provide a formal formulation of the target-driven data refactoring and retrieval. An error-controlled progressive method refactors the original data X = {x1 , . . . , xn } into k progressive segments {s1 , . . . , sk }. Given a user-requested error bound τ during retrieval, the method is supposed to automatically determine the number of fragments required (denoted as nτ ), and use them to reconstruct a decompressed data X ′ = {x′1 , . . . , x′n } which satisfies maxi |xi −P x′i | ≤ τ . The nτ size(si ). total retrieved size under τ is defined as Sτ = i=1 We identify two related but different targets for progressive retrieval based on the diverse needs of scientific applications: (1) error bound and (2) peak Signal-to-Noise ratio (PSNR). An error bound limits the maximal error in the reconstructed data, which is needed in applications that require guaranteed absolute error control, as mentioned in many errorcontrolled lossy compressors [4], [6], [7]. PSNR is defined as 20 log10 (maxi xi − mini xi )/RM SE, where RM SE = pP n ′ 2 i=1 (xi − xi ) /n denotes the root of mean squared errors. It measures the average error over the entire domain, and this is preferred in applications that care about the global error distribution, which is common in fusion energy science and climate studies [29], [30]. While both of them reflect the quality of the reconstructed data, there is no universal solution that optimizes both. Accordingly, we formulate our two targets separately as min Sτ
s.t.,
∥X − X ′ ∥∞ ≤ τ,
for error-bounded retrieval, and max PSNR(X, X ′ )
s.t.,
S = Sτ , ∥X − X ′ ∥∞ ≤ τ,
for PSNR-oriented retrieval. In the rest of the paper, we refer to our designs toward the two targets as the error-bound mode and the PSNR mode, respectively. B. Overview We present an overview of the proposed framework in Figure 1, with blue boxes representing proposed components and cyan boxes indicating optimized components, and gray boxes denoting components adjusted accordingly in scientific data refactoring and progressive retrieval pipelines. Similar to existing methods, our framework has a data refactoring stage, which writes data upon generation, and a retrieval stage, which retrieves data with guaranteed error control for post hoc data analytics. During data refactoring, we leverage and optimize two complementary interpolation schemes and encoding methods to enable an adaptive, flexible refactoring pipeline for diverse targets. We also propose coefficient decomposition, which exploits the commonly overlooked spatial correlation within the decorrelated coefficients to improve retrieval efficiency. During retrieval, we adjust the retrieval
RAW
Adaptive Interpolation
Retrieval Estimation
Proposed components
Coefficient Decomposition
Bitplane Decoding
Coefficient Recomposition
Optimized components
Adaptive Bitplane Encoding
Interpolation
Per-level interpolation
REC
Adjusted components
Fig. 1. Overview of the proposed framework.
estimation, bitplane decoding, and interpolation components to incorporate the adaptive interpolation schemes and the proposed coefficient decomposition. IV. A DAPTIVE I NTERPOLATION We leverage two multilevel interpolation schemes in our framework to adapt to diverse targets, each with distinct characteristics. We refer to them as per-level interpolation and per-region interpolation, respectively, based on how they interpolate data within a level. In this section, we introduce the two schemes and analyze their efficiency with respect to the two targets (error bound and PSNR). Both interpolation schemes decompose the data into levels with different strides in a bottom-up fashion. Levels are used to represent different subsets of data points in an embedded hierarchy, as detailed below. The decomposition starts with a stride of 1 at the finest level (level 0), where all data points are present. For any specific level l, it contains only data points in level l − 1 with even indices along each dimension, i.e., it reduces the resolution by half and doubles the stride in each dimension. The decomposition produces L levels based on user inputs, and terminates at the level L−1, which is typically a coarse representation that accounts for a very small fraction of the data. While the two schemes share the same multilevel decomposition, they feature distinct intra-level interpolation methods for data decorrelation. The procedures for linear interpolation with 3D datasets are illustrated in Figure 2, and they generally extend to any arbitrary dimensionality and higher-order interpolation methods. We describe the two schemes below and highlight our optimizations relative to existing work. Per-level interpolation: The per-level interpolation scheme originates from PMGARD [16], where the intra-level interpolation uses only data points from a higher level. As shown in the top row of Figure 2, it uses the eight data points in the upper level to interpolate all remaining data points in the current level with multilinear interpolation. As such, the data points at the centers of the edges, faces, and cubes are interpolated using the averages of 2, 4, and 8 corresponding upper-level data points, respectively. This scheme has limited error propagation because errors at the current level are completely independent. However, it has relatively low interpolation accuracy because data points at the center of faces and cubes are interpolated using faraway upper-level data points. Compared to PMGARD [16], our per-level interpolation features tricubic interpolation for higher efficiency, and it directly interpolates
z
y
Per-region interpolation
x
Original data points will be processed at level-(l-1). First interpolation direction. Original data points to be processed at level-l. Interpolated by 2 neighboring points. Second interpolation direction. Interpolated by 4 neighboring points. Third interpolation direction. Interpolated by 8 neighboring points. Processed data points after the first direct. interpolation, mapping to Region 3 at level-l. Processed data points after the second direct. interpolation, mapping to Region 2 at level-l. Processed data points after the third direct. interpolation, mapping to Region 1 at level-l.
Fig. 2. Illustration of the decomposition pipelines of per-level linear interpolation and per-region linear interpolation.
the data with strides to avoid the expensive reordering steps for higher throughput. Per-region interpolation: The per-region interpolation scheme is proposed in SZ3 [9] and extended for data decomposition in progressive compression in IPComp [17]. Instead of using fixed data points for interpolation, it adopts a fixed mechanism: each data point in the current level is interpolated using the same formula (e.g., average of two points in the linear example), but with different data. This is typically done by performing separate interpolations along each dimension, as noted in the bottom row of Figure 2. We name all data points interpolated along the same dimension a “region”, so this scheme is called per-region interpolation. Since this scheme uses only adjacent data points for interpolation, it tends to achieve higher accuracy. However, it suffers from error propagation, as errors can propagate across regions at the same level during reconstruction (e.g., from region one to region two when interpolating the centering data point on the top face). Compared to IPComp [17], our per-region interpolation features a bottom-up design and omits the quantization stage, thereby providing more accurate decorrelation using the original data values rather than the quantized ones. We also use a tuner to identify the best-fit interpolation orders across different dimensions, as prior studies [9], [31] show that these orders can affect both interpolation accuracy and throughput. Based on the above analysis, we found that per-level interpolation is preferred for the error-bound mode and per-region interpolation is favorable for the PSNR mode for most cases. As such, we advocate an adaptive interpolation scheme that automatically adjusts based on the target. In practice, we use a tuner to determine which interpolation scheme to use on the fly, as will be detailed in Section VI. V. C OEFFICIENT D ECOMPOSITION Coefficients are the residues between original data values and their interpolated counterparts using the schemes in Section IV. They are generally encoded using bitplane encoding and lossless compression, and their distribution is a key factor in determining the efficiency of progressive compression: a
After Interpolation
Region 1 Region 2 e.g. Region 3
z
y
Region 3
x
Predicted by 2 neighboring points. First prediction direction. Predicted by 4 neighboring points. Second prediction direction. Coefficient points will be processed at coefficient level-(cl-1). Processed coefficient points after the first direct. prediction at coefficient level-cl. Processed coefficient points after the second direct. prediction at coefficient level-cl.
Fig. 3. Illustration of the 2D coefficient decomposition with linear interpolation on 3D data. Coefficients before CoeffDecom MSE = 6.0128𝑒 − 14
Coefficients after CoeffDecom MSE = 4.0978𝑒 − 15
Fig. 4. Slice [:, :, 125] of the 500 × 500 × 250 coefficients in Region 3 at the finest level from the CH4 field of the S3D dataset, decomposed using cubic adaptive interpolation. Left: original coefficients. Right: coefficients after 2D linear CoeffDecom.
large percentage of near-zero coefficients indicates a high percentage of zeros in the most significant bitplanes and thus higher compressibility. Inspired by prior studies [32] that investigated the correlation of quantization indices in error-controlled lossy compressors, we propose an efficient coefficient decomposition method that leverages coefficient correlations for better progressive compression. A. Motivation Following the procedure in [32], we group the coefficients in different regions into separate groups, as shown in Figure 3. We then extract a 2D slice in Region 3 and visualize it in Figure 4 using the CH4 fields of the S3D dataset as an example (see Table III for dataset information). According to this figure, coefficients in this region exhibit high correlation, and similar phenomena are observed in Region 1 and Region 2. In addition, this phenomenon occurs in both per-level interpolation and per-region interpolation, because the residuals tend to be smooth as they are the differences between two smooth sets of values: the original scientific data and the interpolated data. These observations motivate us to further decorrelate the coefficients for better efficiency. B. Methodology Prior studies [32] use a Lorenzo predictor to decorrelate the quantization indices in error-controlled lossy compressors, but this approach cannot generalize to progressive compression. This is because the Lorenzo predictor exhibits unbounded error propagation, which is not a problem for quantization indices with a fixed value (i.e., errors are all 0) but matters
for coefficients whose values vary with the number of retrieved bitplanes. Instead, we propose to further decompose the coefficients using per-level interpolation, as detailed below. The latter half of Figure 3 depicts our coefficient decomposition method on Region 3, with coefficients computed from per-region interpolation, and the same procedure applies for other regions and per-level interpolation. In particular, we treat the 3D coefficients as a stack of multiple 2D slices, and perform 2D per-level multilinear interpolation on each slice. This is done on each region to ensure all the correlations are properly handled. We use 2D interpolation instead of 3D interpolation because such correlations are mainly seen in the hyperplane orthogonal to the interpolation direction, which also aligns with the observations in [32]. We use the per-level scheme and multilinear interpolation because they have minimal error propagation, and we fix the total number of levels to three for the same reason. The decomposed coefficients for the same data are visualized in Figure 4. It is clearly observed that they have more near-zero data points with much smaller mean squared errors. The coefficient decomposition method can be applied at any level, but we observe a noticeable improvement mainly at the finest level. This is because coefficients at the higher levels (1) only take up a small percentage of data (e.g., less than 1/8 for 3D data); and (2) have weaker correlation as they are spatially farther away from each other. As such, we only perform coefficient decomposition at the finest level throughout the paper. C. Coefficient storage and retrieval While coefficient decomposition reduces entropy in the coefficients, it introduces new challenges for retrieval due to interleaved regions and added levels. In what follows, we detail our adjustments to ensure error-controlled retrieval when coefficient decomposition is coupled with the two interpolation schemes in Section IV. a) Per-level interpolation: In the original design, one level of coefficients in per-level interpolation is encoded into bitplanes as a whole. The maximum value of the coefficients is stored as metadata and used to establish 1-to-1 mappings from bitplane to absolute errors, which provide essential information for error-controlled retrieval algorithms (greedy-based in PMGARD [16] and dynamic programming in IPComp [17]). With coefficient decomposition, such mappings cannot be established directly because the actual absolute error at the level is the maximum error across all regions. To address this issue, we couple the greedy-based retrieval algorithm in [16] and a max ordering mechanism to establish such mappings. In particular, we first leverage the greedybased algorithm to determine the loading order of bitplanes and their corresponding error bounds within each region, and then identify a set of similar error bounds across all regions. After that, we apply max ordering to arrange the bitplanes across all the regions at the level. Since the regions contribute equally to the global error bound through the max operation, we merge their orderings by sorting all steps across the three
Region 1
Region 2
Region 1
Region 3
Level 1
Level 1
Level 2
Level 2
Level 3
Level 3
✏1
Region 2
✏2 <latexit sha1_base64="q3bJu/fuXnANxmLI3MeYVZRpBZk=">AAAB8nicbVBNS8NAEN3Ur1q/qh69LBbBU0mKVI9FLx4r2A9IQ9lsJ+3SzW7Y3Qgl9Gd48aCIV3+NN/+NmzYHbX0w8Hhvhpl5YcKZNq777ZQ2Nre2d8q7lb39g8Oj6vFJV8tUUehQyaXqh0QDZwI6hhkO/UQBiUMOvXB6l/u9J1CaSfFoZgkEMRkLFjFKjJX8ASSacSmGjcqwWnPr7gJ4nXgFqaEC7WH1azCSNI1BGMqJ1r7nJibIiDKMcphXBqmGhNApGYNvqSAx6CBbnDzHF1YZ4UgqW8Lghfp7IiOx1rM4tJ0xMRO96uXif56fmugmyJhIUgOCLhdFKcdG4vx/PGIKqOEzSwhVzN6K6YQoQo1NKQ/BW315nXQbda9Zbz5c1Vq3RRxldIbO0SXy0DVqoXvURh1EkUTP6BW9OcZ5cd6dj2VrySlmTtEfOJ8/tSyQ5A==</latexit>
<latexit sha1_base64="s5j3WYozoEfXWJpy7upNTS5zL+8=">AAAB8nicbVBNS8NAEJ3Ur1q/qh69BIvgqSQi1WPRi8cK9gPaUDbbTbt0sxt2J0IJ/RlePCji1V/jzX/jps1Bqw8GHu/NMDMvTAQ36HlfTmltfWNzq7xd2dnd2z+oHh51jEo1ZW2qhNK9kBgmuGRt5ChYL9GMxKFg3XB6m/vdR6YNV/IBZwkLYjKWPOKUoJX6A5YYLpQc+pVhtebVvQXcv8QvSA0KtIbVz8FI0TRmEqkgxvR9L8EgIxo5FWxeGaSGJYROyZj1LZUkZibIFifP3TOrjNxIaVsS3YX6cyIjsTGzOLSdMcGJWfVy8T+vn2J0HWRcJikySZeLolS4qNz8f3fENaMoZpYQqrm91aUToglFm1Iegr/68l/Suaj7jXrj/rLWvCniKMMJnMI5+HAFTbiDFrSBgoIneIFXB51n5815X7aWnGLmGH7B+fgGs6eQ4w==</latexit>
Segment 1
Region 3
Segment 2
Fig. 5. Establishing segment-error bound mappings across three regions using two identified error bounds. A small rectangle represents a bitplane in the corresponding level, and our algorithm merges bitplanes across regions to form segments enforcing a set of target error bounds.
regions in descending order of their cumulative error. When multiple steps from different regions have the same cumulative error, they are merged into a single step. This results in a unified loading order for bitplanes across all regions at the level, along with the corresponding max error after each step. Figure 5 demonstrates how this algorithm works with two sample error bounds ϵ1 and ϵ2 , and it generally extends to any number of error bounds. In particular, our algorithm identifies ϵ1 and ϵ2 as two target error bounds across regions, and then it merges bitplanes across all regions to ensure the merged segment enforces these error bounds. To this end, a list of 1-to1 mappings from segments to error bounds is established and all the error bounds are stored as metadata. During retrieval, such information is used to enable guaranteed error control. This increases the metadata overhead, since one level now records a list of error bounds rather than a single maximum value. Nonetheless, this overhead remains negligible relative to the data size. b) Per-region interpolation: Since errors propagate across regions and levels in the same way, prior works [17] apply the retrieval algorithms at the granularity of regions. This yields 3(L − 1) + 1 regions to process, given L levels in total. Since our regions are the same as those in per-region interpolation, we directly split any region with coefficient decomposition to three new regions, each of which corresponds to a decomposed coefficient level. After that, we can directly apply existing algorithms, such as the one in [17], with only minimal modification to the error estimation mechanisms (i.e., changing the error estimation from cubic interpolation to linear interpolation in the decomposed coefficient regions). VI. I MPLEMENTATION AND O PTIMIZATION In this section, we introduce our detailed implementation of the adaptive progressive compression pipeline. We first present our target-driven tuning for automatic component selection, followed by a detailed algorithm along with tailored performance optimizations. A. Target-driven tuning We leverage a target-driven tuning method to determine the configuration for our data refactoring pipeline online. In particular, we uniformly sample 1% of the original data and evaluate the retrieval efficiency of each option across 9 commonly used tolerances which are 10−3 , 5 × 10−4 ,
10−4 , . . . , 10−7 . To this end, we employ a voting scheme to identify the best-fit configuration, i.e., the one that receives the most votes in the trials. We also conduct an offline study to justify the 1% sampling ratio, comparing the configuration voted from the sample against the one voted from the full data across all fields of all datasets for a single-machine experiment. Experimental results demonstrate that sampled tuning reproduces the full-data decision for 85.0% of the interpolation-scheme selections and 88.1% of the coefficientdecomposition selections, while reducing the tuning overhead from 87.2% of the total refactoring time down to 11.5%. In addition, the residual disagreements are largely benign. Evaluating the selected configurations over a wider range of 16 tolerances, from 10−1 to 10−9 , i.e., beyond the range on which the vote is taken, the resulting bitrate differs by only 1.38% on average over the cases where the two tuners diverge, and by 0.37% when averaged over all cases. At the loose tolerances outside the voting range, neither tuner is optimized, and the sampled tuner occasionally yields a slightly lower bitrate. We therefore adopt 1% sampling for online tuning. TABLE I T UNING AND TUNED CONFIGURATION FOR 3D DATA
Interpolation schemes Coefficient decomposition Encoding
Error-bound mode
PSNR mode
Per-level / Per-region (6) No / 3 directions Generic bitplane
Per-region (6) No / 3 directions Negabinary with XOR
Table I summarizes the tuning configurations in our framework. In particular, we will evaluate different interpolation schemes, with per-region interpolation comprising 6 variants that represent different permutations of the interpolation order. For coefficient decomposition, we evaluate four options, including three for performing decomposition along 3 different directions and one for skipping. The underlined configurations in the table are called tuned, as they consistently exhibit higher efficiency than their alternatives in offline studies. In particular, we found that (1) the generic bitplane encoding is always better than negabinary encoding in the error-bound mode but worse in PSNR mode, and (2) including XOR with negabinary encoding always leads to better efficiency. These findings align with the design choices in prior work [16], [17]. In addition, we found that per-region interpolation is consistently better than per-level interpolation in PSNR mode, and therefore omit the latter for online tuning. B. Algorithm We introduce our refactoring and retrieval algorithm with per-level interpolation for demonstration purposes, and a similar procedure applies to per-region interpolation. Algorithm 1 presents our refactoring algorithm. We initialize the current stride to 1 in the beginning and perform interpolation at the finest level (lines 1-2). Afterward, we interleave the data into different regions to perform coefficient decomposition (lines 3-16). In particular, we first iterate and process coefficients in each region, and we treat them as a stack of multiple slices (lines 5-7). Coefficients in each slice are decomposed into
Algorithm 1 Refactoring with Coefficient Decomposition
Algorithm 2 Retrieval with Coefficient Decomposition
Input: input data X, number of levels L Output: bitplanes {bpi } and metadata {Mi } 1: stride ← 1 /*Initialize stride*/ 2: Π ← CubicInterpolate(X, stride) /*Perform interpolation*/ 3: regions ← Interleave(X − Π, stride) /*Interleave regions*/ 4: bp ← ∅ /*Initialize the set to store region bitplanes*/ 5: for region ∈ regions do 6: /*Iterate and process each region*/ 7: for S ∈ region do 8: /*Iterate and decompose each slice*/ 9: Π1 ←LinearInterpolate(S, 1) /*Interpolate with stride 1*/ 10: Π2 ←LinearInterpolate(S, 2) /*Interpolate with stride 2*/ 11: S1 , S2 , S3 ← ExtractLevel(S − Π1 − Π2 ) /*Extract decomposed coefficients based on their levels*/ 12: R1 ← R1 ∪ S1 , R2 ← R2 ∪ S2 , R3 ← R3 ∪ S3 /*Collect decomposed coefficients based on their levels*/ 13: end for 14: bp ← bp∪ BitplaneEncoding(R1, R2, R3) /*Encode decomposed coefficients on all levels*/ 15: end for 16: bp0 , M0 ←MaxOrdering(bp) /*Establish segment-error bound mapping, see Section V-C*/ 17: for l = 1 → L − 1 do 18: stride ← 2 ∗ stride 19: Π ← CubicInterpolate(X, stride) 20: bpl , Ml ← BitplaneEncoding(X − Π, stride) 21: end for 22: return {bp0 , bp1 , · · · , bpL−1 } ∪ {M0 , M1 , · · · , ML−1 }
Input: requested error bound τ , bitplanes {bpi } and metadata {Mi } Output: reconstructed data X ′
three levels and collected separately (lines 9-12). When all the slices are processed in the region, bitplane encoding is performed on all the levels to produce binary streams (line 14). After all the regions are processed, the max ordering strategy (see Section V-C) is used to merge bitplanes to segments with metadata M0 denoting the segment to error bound mappings (line 16). This completes the coefficient decomposition, and the remaining procedure for refactoring later levels is the same as the one in PMGARD [16], except that cubic interpolation is used to achieve better efficiency. We then present the corresponding retrieval algorithm in Algorithm 2. We start with the largest stride because reconstruction has to start from the coarsest level (line 1). After initializing reconstructed data (line 2), we determine which bitplanes to retrieve using an interpretation algorithm and the metadata. For per-region interpolation, we directly use the dynamic programming-based method in [17], as it delivers better efficiency with negligible performance overhead. For per-level interpolation, we cannot use the method directly because the merged segments have a different format in metadata. As such, we use the best-first search algorithm with pruning and incremental invocation to find the appropriate bitplanes with similar efficiency, but it incurs slightly higher overhead. After the bitplanes are identified and fetched, we reconstruct the data from level L − 1 to 1 as in the existing approach (lines 4-9). At the finest level, we reconstruct the slices within each region one by one and merge them to form the regions (lines 10-20). In the end, we restore the regions to their corresponding locations and add them back to the data after the final interpolation (lines 22-24). We fix the decomposition depth to three levels based on offline studies. Fewer levels leave residual spatial structure
1: stride ← 2L /*Initialize stride*/ 2: X ′ ← 0 /*Initialize reconstructed data*/ 3: {bpi } ← RetrievalInterpreter(τ , {Mi }) /*Determine how many bitplanes to retrieve*/ 4: for l = L − 1 → 1 do 5: stride ← stride/2 6: C ← BitplaneDecoding({bpi }) /*Decode coefficients*/ 7: Π ← CubicInterpolate(X, stride) /*Perform Interpolation*/ 8: X ′ ← Π+ ExpandByStride(C, stride) /*Add coefficients back*/ 9: end for 10: for region ∈ regions do 11: region ← ∅ /*Initialize region*/ 12: R1 , R2 , R3 ← BitplaneDecoding(bp0 ) /*Decode bitplanes for all levels in region*/ 13: for S ∈ region do 14: /*Iterate and recompose each slice*/ 15: S1 , S2 , S3 ←GetSlice(R1 , R2 , R3 ) /*Get slice from decoded region*/ 16: Π1 ←LinearInterpolate(S3 , 1) /*Interpolate level 2*/ 17: Π2 ←LinearInterpolate(S2 , 2) /*Interpolate level 1*/ 18: S ← ExpandByStride(S1 , 4) + ExpandByStride(S2 , 2) + S3 + Π1 + Π2 /*Reconstruct the slice*/ 19: region ← region ∪ S /*Merge slice to region*/ 20: end for 21: end for 22: C ← Reposition(regions) /*Put recomposed coefficient back to correct locations*/ 23: Π ← CubicInterpolate(X ′ , stride) /*Perform interpolation*/ 24: X ′ ← Π + C /*Add recomposed coefficients back*/ 25: return X ′
in the coefficients under-exploited, whereas deeper hierarchies propagate error across more levels. To quantify this trade-off, we compare the three depths on the SCALE dataset at 16 error bounds spanning 10−1 to 10−9 , averaging the bit-rate over 11 fields at each. At every error bound, we take the lowest of the three bit-rates as the reference and measure the excess of each depth over it. For instance, at ϵ = 10−5 in the error-bound mode, three-level yields the lowest bit-rate, while two-level and four-level exceed it by 1.1% and 0.9%, respectively. Table II reports, for each depth, the number of error bounds at which it is the cheapest (#best) together with its average and worst excess across the 16 error bounds. No single depth wins everywhere, but the penalty for a wrong choice is strongly asymmetric. Three-level is the cheapest at 11 of the 16 error bounds in the error-bound mode and at 13 in the PSNR mode, and whenever it loses it does so by at most 0.83%. Four-level gives up as much as 3.71% where it is not preferred, and twolevel never wins at all while costing up to 5.00%. The error bounds favoring four-level also differ between the two modes, so no deeper configuration is uniformly preferable. According to these observations, we adopt three as the decomposition depth throughout the paper. C. Performance optimization While our major goal is to improve the quality of progressive compression, keeping reasonable throughput is equally important, especially for data transfer tasks when refactoring/reconstruction operations lie on the critical paths. This
section details our efforts to optimize the performance of the proposed framework. a) Fastest direction interpolation: Inspired by recent work [31], we found that all interpolations can be enforced to be performed along the fastest-varying direction for better cache efficiency and thus higher performance. In particular, we observe that the per-level linear interpolation described in Section IV is mathematically equivalent to sequentially interpolating the three regions within a level: first predicting round-type points along the x-direction using 2 surrounding triangle vertices, then predicting star-type points along the ydirection using 2 already-updated round points, and finally predicting diamond-type points along the z-direction using 2 already-updated star points (see Figure 2 for reference). This reformulation does not change the interpolation directions or the predicted values; it only changes the order in which data points are visited. Crucially, each region can be traversed along the fastest-varying dimension (i.e., the contiguous dimension in memory), ensuring that all interpolation steps remain cachefriendly without any data reordering. This strategy also applies to per-region interpolation: since each region is interpolated independently, each region naturally corresponds to its own buffer, and the traversal within each region follows the fastest direction as well. b) Optimized generic bitplane encoding: As noted in PMGARD [16], generic bitplane encoding is relatively slow despite its better efficiency in the error-bound mode, because the encoding operation has a branch to check if any bit has been stored for the current data. While this is necessary to ensure correctness and efficiency, it is no longer needed when the sign bits for all data points are recorded. As such, we decompose the algorithm into two sections: the first section retains the original branch when at least one data point has not recorded its sign, and it switches to the second section when all signs are recorded. In the second section, we eliminate such a branch to accelerate the encoding process. This yields a substantial speedup over the original per-bit encoding implementation, since most data chunks complete the sign recordings with only a few bitplanes. We validate the proposed optimization using all 9 fields in the S3D dataset and present the results in Figure 6. It is observed that both optimizations yield a 2× performance improvement over existing implementations for the corresponding operations. VII. E VALUATION We evaluate our method, named ProAICD for Progressive compression with Adaptive Interpolation and Coefficient
Field
Te O2 m p Ve . l. Ve X l. Ve Y l. Z
2
0.00 2O
0
CO
4.16% 0.83% 3.71%
H
1.13% 0.13% 2.06%
4
0 13 3
0.05
CO
5.00% 0.81% 3.63%
PMGARD generic Optimized generic
CH
1.63% 0.14% 0.55%
Throughput (GB/s)
0 11 5
1
Encoding Throughput 0.10
Te O2 m p Ve . l. Ve X l. Ve Y l. Z
Worst
2
Avg.
2O
#best
CO
Worst
H
Avg.
2
4
#best
PMGARD Interp Fastest Direct Interp.
CO
Two-level Three-level Four-level
PSNR mode
Throughput (GB/s)
Error-bound mode
Depth
Decomposition Throughput
CH
TABLE II S ENSITIVITY TO COEFFICIENT DECOMPOSITION DEPTH ON SCALE
Field
Fig. 6. Throughput of decomposition and generic bitplane encoding before and after optimization on the S3D dataset with 9 valid fields
Decomposition, and compare it with three state-of-the-art approaches: PMGARD [16], SZ3-R [28], and IPComp [17] in terms of retrieval efficiency, reconstruction quality, and throughput using five real-world datasets. We also present two typical scientific use cases for end-to-end data transfer and snapshot visualization. A. Experimental Setup 1) Benchmark datasets: We evaluate on five real-world scientific datasets spanning multiple scientific domains: Climate Simulation (CESM [33]), Hydrodynamics Simulation (Miranda [34]), Weather Simulation (SCALE [35]), Combustion Simulation (S3D [36]), and large-scale Turbulence Simulation (JHTDB [37]). The information about these datasets is detailed in Table III. Note that we exclude the Temperature field from SCALE and the N2 and Pressure fields from S3D for the aggregated results, as IPComp encounters numerical overflow on these fields, producing a constant tiny bitrate, no error guarantee, and negative PSNR, and we further report the results on excluded fields separately without IPComp. We only evaluate JHTDB in the parallel data transfer experiments due to its large size. TABLE III DATASETS Dataset
Dimensions
Valid Fields
Type
Size
CESM Miranda SCALE S3D JHTDB
26 × 1800 × 3600 256 × 384 × 384 98 × 1200 × 1200 500 × 500 × 500 4096 × 4096 × 4096
33 7 11 9 1
double double double double double
41.42 GB 1.97 GB 11.57 GB 8.38 GB 512 GB
2) Platform: All experiments are conducted on the Morgan Compute Cluster (MCC) [38], a medium-scale cluster with 100 Gbps InfiniBand HDR interconnect. Each compute node in the system is equipped with 2 AMD EPYC ROME 7702P processors, each with 64 cores and 256 GB of memory. Experiments associated with runtime are evaluated three times, and the average number is reported. 3) Quality assessment: We assess the efficiency of progressive retrieval using the widely adopted rate-distortion graph [39]–[41], which depicts the relationship between bitrate and distortion. Bit-rate represents the average number of bits per data point in the retrieved data and serves as the xaxis, computed as Sτ × 8/n, where Sτ is the retrieved size
j=1 N RM SEj
CH4
CO
CO2
H2O
O2
Temperature
VelocityX
VelocityY
VelocityZ
10−2 10−4 10−6
Error Bound
under requested tolerance τ and n is the total number of data points. We include the metadata size in Sτ for all compressors, since metadata must be loaded to interpret retrieval sizes. To evaluate our two optimization targets, we use two types of rate-distortion plots: bit-rate versus relative error bound, and bit-rate versus PSNR. In the former, curves that are lower and farther left indicate higher efficiency; in the latter, curves that are higher and farther left indicate higher efficiency. Since PMGARD utilizes different encoding methods for the two targets, we use the generic bitplane encoding for PMGARD in the bit-rate versus error bound graphs and negabinary encoding in the bit-rate versus PSNR graphs. SZ3-R iteratively compresses the residuals into 18 snapshots with target error bounds from 10−1 to 10−18 to fully preserve the precision limits of double-precision data. Note that bit-rate and PSNR in Figure 9 and Figure 11 are aggregated using all fields from the same dataset due to the limited space. The aggregated bit-rate is easily computed as the average bit-rate of all fields since they all share the same size in one √dataset. The aggregated PSNR is computed nf , where nf stands for number of as 20 log10 Pnf 2
10−2 10−4 10−6
10−2 10−4 10−6 0
2
4
6
8
0
2
4
6
8
0
2
4
6
8
Bitrate PMGARD
AdatInterp
CoeffDecom
Fig. 7. Ablation study of error-bound mode on the S3D dataset. CH4
CO
CO2
H2O
O2
Temperature
VelocityX
VelocityY
VelocityZ
175 150
RM SE
j fields in one dataset, and N RM SEj = max(Xj )−min((X j )) for j-th field.
125 100 75
B. Ablation Study 175
C. Comparison with State of the Arts We then compare our method with 3 state-of-the-art progressive methods as mentioned earlier. 1) Efficiency: We present our improvement in error-bound mode in Figure 9 and in PSNR mode in Figure 11 using all valid fields of the first four datasets in Table III. Results on the three aforementioned excluded fields are presented in Figure 10 and Figure 12. For error-bound mode, as shown in Figure 9, compression ratios are improved by up to 23.9%, 29.7%, 29.3%, and 42.3% respectively, when compared with the best-performing methods between PMGARD, SZ3-R, and IPComp. ProAICD
150
PSNR
We take the S3D dataset with 9 fields for the ablation study to demonstrate our efficiency gain step by step. For error-bound mode, illustrated in Figure 7, we use PMGARD as the baseline because per-level interpolation originates from PMGARD. Adaptive interpolation (AdatInterp) consistently improves efficiency over PMGARD across all nine fields. Coefficient decomposition (CoeffDecom) further widens this gap, especially for CH4 , CO, CO2 , H2 O, O2 , and Temperature. At worst, it preserves efficiency for VelocityX , VelocityY , and VelocityZ rather than degrading it. As for PSNR mode, as shown in Figure 8, we take IPComp as the baseline since IPComp first applies per-region interpolation into the progressive method. According to Figure 8, we observe that AdatInterp still improves the efficiency by exploring the best-fit interpolation order while IPComp fixes the order, and then CoeffDecom obviously further improves the data quality throughout all nine fields.
125 100 75
175 150 125 100 75 0
2
4
6
8
0
2
4
6
8
0
2
4
6
8
Bitrate IPComp
AdatInterp
CoeffDecom
Fig. 8. Ablation study of PSNR mode on the S3D dataset.
outperforms PMGARD and IPComp across all error bounds on all datasets. SZ3-R achieves a lower bit-rate at certain relatively loose error bounds on the CESM and SCALE datasets; however, this is because those error bounds coincide with the exact target error bounds used in its residual-based compression, giving it a natural advantage at those specific operating points. Even though tighter error bounds such as 10−6 and 10−7 are also directly targeted by SZ3-R, this advantage diminishes as the error bound decreases, due to the inherent redundancy in its residual-based compression scheme. We can also observe similar advantages on excluded fields from Figure 10. For PSNR mode, as illustrated in Figure 11, compression ratios are further optimized by up to approximately 92.5%, 51.1%, 75.6%, and 91.3%, respectively, compared with the best-performing methods. Note that unlike non-progressive methods such as [32], where a given error bound deterministi-
cally produces a specific PSNR, progressive methods retrieve coefficients incrementally, so the same requested error bound does not guarantee the same PSNR across different methods. To enable fair comparison, we apply linear regression between points on each rate-distortion curve and compare bit-rates at matched PSNR values. Our method clearly outperforms all baselines across all datasets (including the excluded fields shown in Figure 12) in terms of efficiency in PSNR mode.
𝐶𝑅 ↑= 23.9% 𝑒𝑏 = 5×10!"
𝐶𝑅 ↑= 92.5% 𝑃𝑆𝑁𝑅 = 58.3
𝐶𝑅 ↑= 51.1% 𝑃𝑆𝑁𝑅 = 77.4
𝐶𝑅 ↑= 75.6% 𝑃𝑆𝑁𝑅 = 64.6
𝐶𝑅 ↑= 91.3% 𝑃𝑆𝑁𝑅 = 118.7
𝐶𝑅 ↑= 29.7% 𝑒𝑏 = 5×10!"
Fig. 11. Baseline comparison of PSNR mode.
𝐶𝑅 ↑= 29.3% 𝑒𝑏 = 5×10!"
PSNR
𝐶𝑅 ↑= 42.3% 𝑒𝑏 = 5×10!"
160 140 120 100 80 60
SCALE: Temperature
0
2
4
6
S3D: N2
175 150 125 100 75
8
PMGARD
0
2
4
Bitrate
SZ3-R
S3D: Pressure
175 150 125 100 75 6
8
0
2
4
6
8
ProAICD
Fig. 9. Baseline comparison of error-bound mode. Fig. 12. Baseline comparison of PSNR mode on excluded fields.
SCALE: Temperature
S3D: N2
S3D: Pressure
Error Bound
10 2 10 4 10 6 0
2
4
6
8 0
PMGARD
2
4
Bitrate
SZ3-R
6
8 0
2
4
6
8
ProAICD
Fig. 10. Baseline comparison of error-bound mode on excluded fields.
2) Performance: We report the average refactor time and average reconstruction time of the four methods in Table IV and Table V, using three requested relative error bounds of 10−2 , 10−4 , and 10−6 as an example for reconstruction. For refactoring, PMGARD’s PSNR mode is the fastest across all four datasets, followed by IPComp. Our PSNR mode achieves comparable refactoring time to IPComp, while our error-bound mode is moderately slower due to the additional overhead of tuning, coefficient decomposition, and ordering. SZ3-R is the slowest by a significant margin, taking 4×–10× longer than the other methods, as it iteratively compresses the residuals at each target error bound from 10−1 down to 10−18 . For reconstruction, IPComp is the fastest at τ = 10−4 and 10−6 on CESM, SCALE, and S3D. Note that the reported time corresponds to a single retrieval request at the given tolerance rather than an accumulation over progressively refined requests. Across all methods, reconstruction time increases as the tolerance tightens, since more data must be retrieved
and decompressed; the increase is sharpest for SZ3-R, whose representation is a chain of residual snapshots, so a request at τ must retrieve and decompress every snapshot down to τ , accumulating cost with the length of the chain rather than with the requested precision alone. This trend holds across all three tolerances and all four datasets, confirming that it is not specific to a single operating point. Separately, both IPComp and SZ3-R exhibit notably slow reconstruction on the Miranda dataset, despite it being the smallest dataset, and this ranking also holds across all three tolerances. We observe a similar pattern in both single-machine and parallel experiments on the JHTDB dataset, suggesting that these two methods may not perform well on turbulence simulation datasets in general. In terms of bit-rate, the comparison against SZ3-R is governed by the same accumulation. At the loose end, τ = 10−2 , only the first two components are needed, and it attains the lowest bit-rate on all four datasets, as a short residual chain is hard to beat when only a coarse approximation is requested. At τ = 10−4 , its advantage is already confined to CESM and SCALE, while our error-bound mode is the lowest on Miranda and S3D. Once the chain lengthens further, the accumulated cost of storing successive residuals outweighs the advantage entirely: at τ = 10−6 our approach attains the lowest bitrate on all four datasets, whereas SZ3-R becomes the most expensive, e.g., BR = 9.12 versus BR = 6.42 on CESM
TABLE IV AVERAGE REFACTOR TIME ( IN SECONDS ) OF DIFFERENT PROGRESSIVE APPROACHES USING ALL FIELDS FROM THE SAME DATASET
Datasets Method CESM Miranda SCALE EB 19.20 7.14 21.86 ProAICD PSNR 17.05 3.96 14.63 EB 32.41 9.90 27.24 PMGARD PSNR 12.36 2.82 10.36 SZ3-R – 106.59 42.51 93.18 IPComp – 14.77 4.77 12.17
S3D 18.37 14.49 33.54 9.91 81.48 11.39
TABLE V AVERAGE RECONSTRUCTION TIME ( IN SECONDS ) AND AVERAGE BIT- RATE (BR) OF DIFFERENT PROGRESSIVE APPROACHES USING ALL FIELDS FROM THE SAME DATASET UNDER TOLERANCES τ = 10−2 , 10−4 , AND 10−6 Method
ProAICD PMGARD SZ3-R IPComp ProAICD PMGARD SZ3-R IPComp ProAICD PMGARD SZ3-R IPComp
EB PSNR EB PSNR – – EB PSNR EB PSNR – – EB PSNR EB PSNR – –
CESM Miranda SCALE S3D Time BR Time BR Time BR Time BR τ = 10−2 4.09 0.25 0.98 0.14 3.32 0.15 3.43 0.05 2.82 0.36 0.75 0.21 2.35 0.32 2.59 0.07 4.04 0.34 0.78 0.20 2.96 0.16 1.80 0.10 3.17 0.44 0.68 0.32 2.41 0.28 1.95 0.15 2.11 0.16 2.48 0.07 1.71 0.04 1.44 0.02 3.13 0.37 5.06 0.21 2.92 0.30 2.34 0.08 τ = 10−4 7.12 2.14 1.51 1.17 5.98 2.17 4.82 0.48 3.55 2.36 0.92 1.26 2.99 2.48 3.06 0.56 8.04 2.86 1.53 1.92 6.15 2.56 4.25 1.65 4.33 3.05 0.91 2.06 3.42 2.77 2.83 1.87 4.95 1.92 4.74 1.48 4.12 1.69 3.06 0.88 3.14 2.45 5.06 1.36 2.89 2.51 2.27 0.66 τ = 10−6 12.21 6.42 2.28 3.28 10.25 6.43 7.22 2.50 4.39 6.42 1.10 3.28 3.65 6.43 3.69 2.54 13.73 7.35 2.44 4.56 10.92 7.05 8.39 5.75 5.37 7.53 1.12 4.71 4.30 7.17 3.58 5.91 11.29 9.12 8.46 6.73 9.02 9.02 7.36 7.48 3.69 6.54 5.43 3.49 3.29 6.73 2.74 3.02
D. Use case for eb mode: end-to-end data transfer We demonstrate how our method can benefit scientific data management using a practical use case of remote data transfer. Here is the experiment setup: suppose the refactored data is stored at the Frontier Supercomputer at Oak Ridge Leadership Computing Facilities [1], and the request for data is initiated from MCC [38] with a frequently used tolerance τ = 10−4 . The data is transferred via Globus [42], the leading research cyberinfrastructure widely used for scientific data sharing and management. This setup reflects a common scenario in which domain scientists must download data from remote servers to local machines for post hoc data analysis. To illustrate how the progressive methods perform with increasing data sizes, we conduct a weak-scaling experiment using 128, 256, 512, and 1024 cores, respectively. In all the experiments, each core
End-to-End Data Transfer Time
100
Total Time (s)
compared with ProAICD. Overall, our error-bound mode achieves the lowest bit-rate on the majority of datasets at the cost of slightly longer refactoring and reconstruction time, while our PSNR mode maintains competitive speed with consistently lower bit-rate than all three baselines.
ProAICD Retrieval
SZ3-R Retrieval
ProAICD Transfer
SZ3-R Transfer
PMGARD Retrieval
IPComp Retrieval
PMGARD Transfer
IPComp Transfer
80 60 40 20 0
128 (64GB)
256 (128GB)
512 (256GB)
1024 (512GB)
Number of Cores (Data Size)
Fig. 13. End-to-end data transfer time using JHTDB dataset. Transferring the original data of 128, 256, 512, and 1024 cores takes 189, 373, 735, and 1342 seconds, respectively.
is processing a fixed size (256 × 512 × 512) of data from the JHTDB dataset, so the total size of the processed data increases with the number of cores. We plot the end-to-end data transfer time (which includes both retrieval time and transfer time) in Figure 13. During retrieval, SZ3-R yields the smallest retrieved size because τ = 10−4 is one of its targeted error bounds. However, this advantage is offset by the longest retrieval time among all methods as previously observed on the Miranda dataset, both SZ3-R and IPComp exhibit slow reconstruction on turbulence simulation datasets, and JHTDB falls into the same category. Moreover, IPComp yields the second largest retrieved size, only smaller than PMGARD, further limiting its transferring performance. Our method yields the second smallest retrieved size, only slightly larger than SZ3-R, mainly because it leverages the most correlation within the dataset using adaptive interpolation and coefficient decomposition techniques. Meanwhile, our method achieves the second shortest retrieval time, only slightly slower than PMGARD. As a result, our method takes the shortest end-to-end data transfer time in experiments with 256, 512, and 1024 cores. Specifically, performance gain is up to 1.26× over the best existing method and 18.17× over vanilla transfer of the original data. E. Use case for PSNR mode: data visualization We illustrate the practical impact of progressive retrieval on scientific data visualization in Figure 14, using the Temperature field from the CESM dataset as an example. For the remaining progressive methods, we cap the retrieval bitrate at approximately 0.5 to ensure a fair comparison. Notably, our method produces the most faithful reconstruction at a bitrate of only 0.37 and is visually nearly indistinguishable from the original, whereas the other methods show noticeable artifacts or over-smoothing even at higher bitrates. In other words, our method retrieves substantially less data from disk yet yields
better visual quality, a property that is especially valuable in interactive exploration of large-scale scientific datasets. Full Scale
Region Scale
Original
Science Foundation (NSF) under Grant OAC-2628470, OAC2628471, OAC-2628472, OAC-2311757, and OAC-2144403. This research used resources of the Oak Ridge Leadership Computing Facility (OLCF), which is a DOE Office of Science User Facility. R EFERENCES
IPComp, 𝐵𝑅 = 0.49 𝑃𝑆𝑁𝑅 = 72.5 𝑅𝑂𝐼 𝑃𝑆𝑁𝑅 = 32.1
PMGARD, 𝐵𝑅 = 0.51 𝑃𝑆𝑁𝑅 = 66.6 𝑅𝑂𝐼 𝑃𝑆𝑁𝑅 = 26.3
SZ3-R, 𝐵𝑅 = 0.80 𝑃𝑆𝑁𝑅 = 66.8 𝑅𝑂𝐼 𝑃𝑆𝑁𝑅 = 28.2
ProAICD, 𝐵𝑅 = 0.37 𝑃𝑆𝑁𝑅 = 87.2 𝑅𝑂𝐼 𝑃𝑆𝑁𝑅 = 49.6
Fig. 14. Visual comparison of reconstructed data from different progressive approaches on the CESM dataset, field Temperature, at the slice [11, 1000:1200, 1480:1680]. Each panel shows the method name, achieved bitrate, global PSNR over the full field, and ROI PSNR computed within the displayed region. The left and right colorbars represent the color scales in the full slice and the displayed region, respectively.
VIII. C ONCLUSION In this paper, we design an adaptive progressive compression framework to enable efficient retrieval of scientific data toward diverse targets. Our approach improves both retrieval efficiency and reconstruction quality by introducing a highly adaptive design with respect to interpolation schemes and encoding methods, and by employing a novel coefficient decomposition method that exploits the commonly overlooked correlation among decorrelated coefficients. Experimental evaluations demonstrate that the error-bound mode of our method delivers up to 42.3% improvement in compression ratio with comparable throughput over state-ofthe-art approaches under the same requested tolerance, leading to up to 1.26× speedup in end-to-end data transfer. The PSNR mode further achieves up to 92.5% improvement under the same PSNR with comparable throughput, and produces the best visualization quality with the least retrieved data. In future work, we plan to explore novel algorithms to further exploit correlations within both raw data and decorrelated coefficients, and to extend our framework to GPUs to further improve throughput. ACKNOWLEDGMENT The research is supported in part by the U.S. Department of Energy (DOE) RAPIDS-3 SciDAC and Sirius-2 projects under contract number DE-AC05-00OR22725, and National
[1] “Frontier exascale supercomputer,” https://www.olcf.ornl.gov/frontier. [2] P. Yeung, K. Ravikumar, R. Uma-Vaideswaran, D. L. Dotson, K. R. Sreenivasan, S. B. Pope, C. Meneveau, and S. Nichols, “Small-scale properties from exascale computations of turbulence on a 32 7683 periodic cube,” Journal of Fluid Mechanics, vol. 1019, p. R2, 2025. [Online]. Available: https://doi.org/10.1017/jfm.2025.10493 [3] S. Lakshminarasimhan, N. Shah, S. Ethier, S.-H. Ku, C.-S. Chang, S. Klasky, R. Latham, R. Ross, and N. F. Samatova, “Isabela for effective in situ compression of scientific data,” Concurrency and Computation: Practice and Experience, vol. 25, no. 4, pp. 524–540, 2013. [Online]. Available: https://doi.org/10.1002/cpe.2887 [4] D. Tao, S. Di, Z. Chen, and F. Cappello, “Significantly improving lossy compression for scientific data sets based on multidimensional prediction and error-controlled quantization,” in 2017 IEEE International Parallel and Distributed Processing Symposium. IEEE, 2017, pp. 1129–1139. [Online]. Available: https://doi.org/10.1109/IPDPS.2017.115 [5] P. Lindstrom and M. Isenburg, “Fast and efficient compression of floating-point data,” IEEE transactions on visualization and computer graphics, vol. 12, no. 5, pp. 1245–1250, 2006. [Online]. Available: https://doi.org/10.1109/TVCG.2006.143 [6] 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 [7] M. Ainsworth, O. Tugluk, B. Whitney, and S. Klasky, “Multilevel techniques for compression and reduction of scientific data—the univariate case,” Computing and Visualization in Science, vol. 19, no. 5-6, pp. 65–76, 2018. [Online]. Available: https://doi.org/10.1007/ s00791-018-00303-9 [8] X. Liang, S. Di, D. Tao, S. Li, S. Li, H. Guo, Z. Chen, and F. Cappello, “Error-controlled lossy compression optimized for high compression ratios of scientific datasets,” in 2018 IEEE International Conference on Big Data. IEEE, 2018, pp. 438–447. [Online]. Available: https://doi.org/10.1109/BigData.2018.8622520 [9] 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 2021 IEEE 37th International Conference on Data Engineering (ICDE). IEEE, 2021, pp. 1643–1654. [Online]. Available: https://doi.org/10.1109/ICDE51399.2021.00145 [10] 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, 2022. [Online]. Available: https://doi.org/10.1109/TBDATA.2022.3201176 [11] Ainsworth, Mark and Tugluk, Ozan and Whitney, Ben and Klasky, Scott, “Multilevel techniques for compression and reduction of scientific data—the multivariate case,” SIAM Journal on Scientific Computing, vol. 41, no. 2, pp. A1278–A1303, 2019. [Online]. Available: https://epubs.siam.org/doi/10.1137/18M1166651 [12] Ainsworth, Mark and Tugluk, Ozan and Whitney, Ben and Klasky, Scott, “Multilevel techniques for compression and reduction of scientific data-quantitative control of accuracy in derived quantities,” SIAM Journal on Scientific Computing, vol. 41, no. 4, pp. A2146–A2171, 2019. [Online]. Available: https://doi.org/10.1137/18M1208885 [13] X. Liang, B. Whitney, J. Chen, L. Wan, Q. Liu, D. Tao, J. Kress, D. R. Pugmire, M. Wolf, N. Podhorszki, and S. Klasky, “Mgard+: Optimizing multilevel methods for error-bounded scientific data reduction,” IEEE Transactions on Computers, 2021. [Online]. Available: https://doi.org/10.1109/TC.2021.3092201 [14] T. H. Group, “The hdf5 library & file format,” https://www.hdfgroup. org/solutions/hdf5/, online. [Online]. Available: https://www.hdfgroup. org/solutions/hdf5/
[15] W. F. Godoy, N. Podhorszki, R. Wang, C. Atkins, G. Eisenhauer, J. Gu, P. Davis, J. Choi, K. Germaschewski, K. Huck et al., “Adios 2: The adaptable input output system. a framework for high-performance data management,” SoftwareX, vol. 12, p. 100561, 2020. [Online]. Available: https://doi.org/10.1016/j.softx.2020.100561 [16] X. Liang, Q. Gong, J. Chen, B. Whitney, L. Wan, Q. Liu, D. Pugmire, R. Archibald, N. Podhorszki, and S. Klasky, “Errorcontrolled, progressive, and adaptable retrieval of scientific data with multilevel decomposition,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. ACM, 2021, pp. 1–13. [Online]. Available: https: //doi.org/10.1145/3458817.3476179 [17] Z. Yang, S. Di, L. Zhang, R. Li, X. Li, J. Huang, J. Liu, F. Cappello, and K. Zhao, “Ipcomp: Interpolation based progressive lossy compression for scientific applications,” in Proceedings of the 34th International Symposium on High-Performance Parallel and Distributed Computing. ACM, 2025. [Online]. Available: https://dl.acm.org/doi/10.1145/3731545.3731578 [18] H. Bhatia, D. Hoang, N. Morrical, V. Pascucci, P.-T. Bremer, and P. Lindstrom, “Amm: Adaptive multilinear meshes,” IEEE Transactions on Visualization and Computer Graphics, vol. 28, no. 6, pp. 2350–2363, 2022. [Online]. Available: https://doi.org/10.1109/TVCG.2022.3165392 [19] G. K. Wallace, “The jpeg still picture compression standard,” IEEE transactions on consumer electronics, vol. 38, no. 1, pp. xviii–xxxiv, 1992. [Online]. Available: https://doi.org/10.1109/30.125072 [20] C. Christopoulos, A. Skodras, and T. Ebrahimi, “The jpeg2000 still image coding system: an overview,” IEEE transactions on consumer electronics, vol. 46, no. 4, pp. 1103–1127, 2000. [Online]. Available: https://doi.org/10.1109/30.920468 [21] J. P. Clyne, E. Bethel, H. Childs, and C. Hansen, “Progressive data access for regular grids.” 2012. [Online]. Available: https: //doi.org/10.1201/b12985-11 [22] W. Li, Q. Gong, X. Wu, J. Chen, Q. Liu, X. He, N. Podhorszki, S. Klasky, and X. Liang, “Qpror: An efficient framework for quantity-of-interest based progressive retrieval with guaranteed error control,” in Proceedings of the 35th International Symposium on High-Performance Parallel and Distributed Computing, 2026. [Online]. Available: https://doi.org/10.1145/3806645.3807579 [23] P. Deutsch, “Gzip file format specification version 4.3,” 1996. [Online]. Available: https://doi.org/10.17487/RFC1952 [24] Y. Collet, “Zstandard - real-time data compression algorithm,” http: //facebook.github.io/zstd/, online. [Online]. Available: http://facebook. github.io/zstd/ [25] F. Alted, “Blosc compressor,” http://blosc.org/, online. [Online]. Available: http://blosc.org/ [26] 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 [27] 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 [28] V. A. Magri and P. Lindstrom, “A general framework for progressive data compression and retrieval,” IEEE Transactions on Visualization and Computer Graphics, 2023. [Online]. Available: https://doi.org/10. 1109/TVCG.2023.3327186 [29] Q. Gong, X. Liang, B. Whitney, J. Y. Choi, J. Chen, L. Wan, S. Ethier, S.-H. Ku, R. M. Churchill, C.-S. Chang, M. Ainsworth, O. Tugluk, T. Munson, D. Pugmire, R. Archibald, and S. Klasky, “Maintaining trust in reduction: preserving the accuracy of quantities of interest for lossy compression,” in Smoky Mountains Computational Sciences and Engineering Conference. Springer, 2021. [Online]. Available: https://doi.org/10.1007/978-3-030-96498-6 2 [30] A. H. Baker, D. M. Hammerling, S. A. Mickelson, H. Xu, M. B. Stolpe, P. Naveau, B. Sanderson, I. Ebert-Uphoff, S. Samarasinghe, F. D. Simone, F. Carbone, C. N. Gencarelli, J. M. Dennis, J. E. Kay, and P. Lindstrom, “Evaluating lossy data compression on climate simulation data within a large ensemble,” Geoscientific Model Development, vol. 9, no. 12, pp. 4381–4403, 2016. [Online]. Available: https://doi.org/10.5194/gmd-9-4381-2016 [31] J. Liu, S. Di, K. Zhao, X. Liang, S. Jin, Z. Jian, J. Huang, S. Wu, Z. Chen, and F. Cappello, “High-performance effective scientific error-bounded lossy compression with auto-tuned multi-
component interpolation,” Proceedings of the ACM on Management of Data, vol. 2, no. 1, pp. 1–27, 2024. [Online]. Available: https://doi.org/10.1145/3639259 [32] P. Jiao, S. Di, M. Xia, X. Wu, J. Liu, X. Liang, and F. Cappello, “Improving the efficiency of interpolation-based scientific data compressors with adaptive quantization index prediction,” in 2025 IEEE International Parallel and Distributed Processing Symposium (IPDPS), 2025, pp. 974–986. [Online]. Available: https://doi.org/10. 1109/IPDPS64566.2025.00091 [33] J. E. Kay, C. Deser, A. Phillips, A. Mai, C. Hannay, G. Strand, J. M. Arblaster, S. Bates, G. Danabasoglu, J. Edwards, M. Holland, P. Kushner, J.-F. Lamarque, D. Lawrence, K. Lindsay, A. Middleton, E. Munoz, R. Neale, K. Oleson, L. Polvani, and M. Vertenstein, “The community earth system model (cesm) large ensemble project: A community resource for studying climate change in the presence of internal climate variability,” Bulletin of the American Meteorological Society, vol. 96, no. 8, pp. 1333–1349, 2015. [Online]. Available: https://doi.org/10.1175/BAMS-D-13-00255.1 [34] “Miranda application,” https://wci.llnl.gov/simulation/computer-codes/ miranda. [35] G.-Y. Lien, T. Miyoshi, S. Nishizawa, R. Yoshida, H. Yashiro, S. A. Adachi, T. Yamaura, and H. Tomita, “The near-real-time scale-letkf system: A case of the september 2015 kanto-tohoku heavy rainfall,” SOLA, vol. 13, pp. 1–6, 2017. [Online]. Available: https://doi.org/10.2151/sola.2017-001 [36] J. Chen, “S3d-legion: An exascale software for direct numerical simulation of turbulent combustion with complex multicomponent chemistry,” in Exascale Scientific Applications. Chapman and Hall/CRC, 2017, pp. 257–278. [Online]. Available: https://doi.org/10. 1201/b21930-12 [37] D. Rosenberg, A. Pouquet, R. Marino, and P. D. Mininni, “Evidence for bolgiano-obukhov scaling in rotating stratified turbulence using high-resolution direct numerical simulations,” Physics of Fluids, vol. 27, no. 5, p. 055105, 05 2015. [Online]. Available: https://doi.org/10.1063/1.4921076 [38] “Morgan Compute Cluster,” https://docs.ccs.uky.edu. [39] D. Tao, S. Di, X. Liang, Z. Chen, and F. Cappello, “Optimizing lossy compression rate-distortion from automatic online selection between sz and zfp,” IEEE Transactions on Parallel and Distributed Systems, vol. 30, no. 8, pp. 1857–1871, 2019. [Online]. Available: https://doi.org/10.1109/TPDS.2019.2894404 [40] P. Jiao, S. Di, H. Guo, K. Zhao, J. Tian, D. Tao, X. Liang, and F. Cappello, “Toward quantity-of-interest preserving lossy compression for scientific data,” Proceedings of the VLDB Endowment, vol. 16, no. 4, pp. 697–710, 2022. [Online]. Available: https://doi.org/10.14778/ 3574245.3574255 [41] X. Wu, Q. Gong, J. Chen, Q. Liu, N. Podhorszki, X. Liang, and S. Klasky, “Error-controlled progressive retrieval of scientific data under derivable quantities of interest,” in SC24: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2024, pp. 1–16. [Online]. Available: https: //doi.org/10.1109/SC41406.2024.00092 [42] I. Foster and C. Kesselman, “Globus: A metacomputing infrastructure toolkit,” The International Journal of Supercomputer Applications and High Performance Computing, vol. 11, no. 2, pp. 115–128, 1997. [Online]. Available: https://doi.org/10.1177/109434209701100205