Error-controlled Progressive Retrieval of Scientific Data under Derivable Quantities of Interest

Our progressive retrieval framework is designed to extract only the “necessary” amount of progressive fragments and guarantee that the user-prescribed error tolerance in QoIs derived from the reconstructed data is met. The capability of estimating the errors in QoIs is crucial for the trustability of the reconstructed data, and it can accelerate the process of reading data from low-tier storage or remote central repository by minimizing the data volume. Below, we define the requested functionalities in progressive compressors and the derivable QoIs targeted in this paper.

Despite the fact that there are only seven families of QoIs listed above, their compositions cover a broad range of functions commonly used in scientific studies. Below, we showcase how to break down the user-requested QoIs into the basis of derivable QoIs using the non-proprietary data generated by a Computational Fluid Dynamics (CFD) code from GE. The CFD simulation produces velocity $V_{x}$, $V_{y}$, $V_{z}$, pressure $P$, and density $D$ on a set of unstructured meshes. We linearize the data into a one-dimensional array, choosing six QoIs used in the posthoc analyses with the detailed equations listed below, as they will be used in our experimental evaluation later on:

Here, $R=287.1$, $\gamma=1.4$, $mi=3.5$, $\mu_{r}=1.716e-5$, $T_{r}=273.15$, and $S=110.4$ are different constant values, and the input and intermediate variables are bolded. Using $PT$ defined in Equation (5) as an example, the formula can be decomposed into the multiplication of $P$ and $(1+\frac{\gamma}{2}*Mach*Mach)^{mi}$, where the latter is a composition of the square root function and a polynomial of $Mach$.

We illustrate the workflow of the proposed QoI-preserving progressive data retrieval framework in Fig. 1. Our key contribution lies in the retrieval procedure, which iteratively refines the reduced approximation till the estimate errors in QoIs drop below user-prescribed bounds. Specifically, we develop (1) a QoI error estimator that provides an upper bound of the errors in QoI given the reconstructed data and its $L^{\infty}$ error bounds, and (2) a primary data (PD) error assigner which estimates and prescribes the bounds used for the next round of data retrieval. The progressive compression procedure can be performed using the existing error-controlled progressive compression frameworks [4, 15, 16], where data is refactored and compressed into multi-precision fragments.

The overall data retrieval pipeline can be summarized as follows. Firstly, an analytic task requests a set of QoIs and the desired error tolerance. This request is processed by the PD error-bound assigner (module 1), which gauges the error bounds on each primary data field used by the first round of retrieval. Such error bounds will be sent to a progressive retriever (module 2), which extracts progressive segments from the most to the least significant until the errors in the errors in the reconstructed data reach below the requested bounds. Data will be incrementally recomposed using the newly arrived progressive segments and then fed into the QoI error estimator, along with the error bounds used during retrieval, to estimate the upper bounds of QoI errors under the current data representations (module 3). If the estimated QoI errors are less than the requested tolerances, the data is provisioned for the analyses; otherwise, the current data representations, along with the derived QoI errors, will be forwarded to the PD error bound assigner to estimate the error bounds used in the next round of data retrieval. The pipeline repeats these steps till the targeted QoI error bounds are reached, or a full-fidelity data representation is retrieved. Due to their progressive nature, only incremental portions of the data need to be retrieved in the later requests, which promises high efficiency in managing the data movement from storage systems to applications.

We leverage the widely used rate-distortion curves [15, 41, 42] to evaluate the efficiency of our approach. The X-axis in the curve is bitrate, which represents the average number of bits in the compressed format. It is analogous to the compression ratio in single-snapshot compression, and can be computed by the retrieved data size divided by the number of elements. We use relative QoI errors as our distortion metric for the Y-axis, which is computed by the maximal absolute error of QoI divided by its respective value range. An easy way to compare multiple rate-distortion curves is to fix either the X value or Y value: in the former case, one can compare the errors of different approaches based on the same retrieved data size; in the latter case, one can compare the size of retrieved data under the same quality.

Here, we assume the original value is $x^{\prime}$ will satisfy the error bound constraint $\lvert x^{\prime}-x\rvert\leq\epsilon$, then we have $\lvert f(x^{*})-f(x)\rvert_{\infty}\leq\sup_{\lvert x^{\prime}-x\rvert\leq\epsilon}{\lvert f(x^{\prime})-f(x)\rvert}_{\infty}=\Delta(f,x,\epsilon)$. Note that $\Delta(f,x,\epsilon)$ only relies on the reconstructed data and the error bound used during data retrieval. Below, we present the theorems used for estimating $\Delta(f,x,\epsilon)$ given several univariate QoI functions.

Note that Theorem 3 does not apply to the case of $\epsilon>\lvert x+c\rvert$, as it may lead to an infinitesimal value of $\lvert x+\xi+c\rvert$, causing the errors in QoI unable to be bounded. Such a case can be avoided by only choosing $\epsilon<\lvert x+c\rvert$ during data retrieval.

Assume $\xi_{i}=x_{i}^{\prime}-x_{i}$ and $x^{\prime}_{i}\in[x_{i}-\epsilon_{i},x_{i}+\epsilon_{i}]$ following the error bound constraints, then we have the following theorems.

This subsection applies the theories of addition, scalar multiplication, and composition for the QoI functions presented in the previous subsection to broaden the range of QoIs that can be preserved during retrieval. Specifically, we have the following theorems for composite QoIs.

Although the above proofs have been conducted on univariate QoIs, the same theorems apply to the additive, multiplicative, and composite operations in multivariate QoIs. Additionally, we derive the error estimators for the composition of univariate QoIs and multivariate QoIs and summarize them in the two lemmas below. We omit the details of the proofs due to limited space. The proofs are similar to the procedure in Theorem 9.

These composite theorems and lemmas greatly extend our flexibility, allowing for progressively retrieving and bounding errors in a variety of QoIs. For instance, multiplications of multiple variables in the form of $\Pi{x_{i}}$ can be preserved by iteratively leveraging the multiplication theory (Theorem 5) and composite property (Theorem 9). Errors in a general polynomial in the form of $\sum{a_{i}x^{i}}$ can be upper bounded by applying the additive (Theorem 7), multiplicative (Theorem 8) properties, and the error estimation of power functions (Theorem 1).

Here, we showcase how to estimate the $V_{\textrm{total}}$ in the GE case study using the proposed methods. Let $x_{1},x_{2},x_{3}$ denote $V_{x}$, $V_{y}$, $V_{z}$. The $V_{\textrm{total}}$, as denoted in Equation (1), can be formulated as the composition of a univariate function $f_{1}(x)=\sqrt{x}$ with the composition of a multivariate function $g_{1}(x_{1},x_{2},x_{3})=x_{1}+x_{2}+x_{3}$ and an univariate function $f_{2}(x)=x^{2}$, which yields $V_{total}=f_{1}(g_{1}(f_{2}(x_{1}),f_{2}(x_{2}),f_{2}(x_{3})))$. We estimate the upper bound of errors in $f_{2}$, $g_{1}$, and $f_{1}$ sequentially. First, the upper bound of errors in $f_{2}(x_{i})$ ($i=1,2,3$) can be estimated using Theorem 1. This forms the error bound vector $\bm{\epsilon_{f2}}=(\Delta(f_{2},x_{1},\epsilon_{1}),\Delta(f_{2},x_{2},\epsilon_{2}),\Delta(f_{2},x_{3},\epsilon_{3}))^{T}$. Meanwhile, the value of $f_{2}(x_{i})$ can be computed to obtain the new value vector $\bm{x_{f2}}=(f_{2}(x_{1}),f_{2}(x_{2}),f_{2}(x_{3}))$. After that, $\bm{\epsilon_{f2}}$ and $\bm{x_{f2}}$ will be used to compute $\epsilon_{g1}=\Delta(g_{1}\circ f_{2},\bm{x_{f2}},\bm{\epsilon_{f2}})$ using Theorem 4. At last, we will compute $x_{g1}=g_{1}(f_{2}(x_{1}),f_{2}(x_{2}),f_{2}(x_{3}))$ to derive the final QoI error bound $\Delta(f_{1}\circ g_{1}\circ f_{2},x_{g1},\epsilon_{g1})$ using Theorem 2.

While using QoIs in GE as a demonstrative example, our theories are extendable to other scientific applications due to the following reasons. First, some common QoIs in the paper can be directly used by other applications (e.g., total velocity in climatology and cosmology). Second, the set of basis QoIs can composite diverse and complex functions such as multivariate polynomials and rational functions, which cover a broad range of QoIs, including molar concentration multiplications in combustion. Third, our theory can extend to new operators with derivable error control (e.g., isosurface [21]). This demonstrates the genericity of the proposed QoI-preserving theory.

Our pipeline consists of two stages: a data refactoring stage, which transforms the original data into multi-precision segments for storage, and a data retrieval stage, which fetches and recomposes data till it reaches user-specified QoI tolerances. We omit the discussion on progressive refactoring as it’s a direct application of the existing techniques. Generally speaking, for every data field, the refactoring stage produces a set of multi-precision segments and the corresponding metadata.

The proposed QoI-preserved data retrieval pipeline is presented in Algorithm 2, assuming that the retrieval starts from scratch. Lines 1-6 show the initialization of the progressive data representation and error bounds used in the first round of retrieval. The while loop starting from line 7 presents the iterative procedure used for QoI-preserving data retrieval that the data representations are gradually refined and then checked for QoI errors till all QoI tolerances are satisfied or a full-fidelity data representation has been retrieved. The progressive_construct function takes the newly retrieved multi-precision segments and recomposes the current data representation to a more accurate approximation. Lines 13-23 estimate the errors of QoIs derived from the reconstructed data using the theorems and lemmas discussed in Section IV, and compare them against user-requested error tolerance. When the requested tolerances are not met, we record the maximal estimated errors as well as their corresponding locations in the data space. Such information will be used for optimizing the error bound assigned for the next round of data retrieval (lines 17 - 22).

Regarding the PD error bound used for data retrieval, the initialization and iterative refinement stages adopt separate error assignment algorithms. At the initialization (line 5) stage, we adopt the Algorithm 3: when a data field is utilized by multiple QoIs, its error bound will be determined by the minimal relative tolerance among all the requested QoIs involving this field. At the iterative refinement (line 26) stage, we adopt a uniform error-tightening strategy detailed in Algorithm 4, and prioritize the evaluation of the data point that generates the largest QoI errors in the previous evaluation. If the estimated QoI errors exceed the tolerance, we reduce the error bounds of all the variables involved in the computation for this QoI by a constant factor $c$ ($c=1.5$ in our implementation), iteratively retrieve the data, then estimate the QoI errors under the new error bounds and the reconstruction. Note that we execute the QoI error estimation only on the data point with the largest QoI error at the iterative refinement stage, which decreases the number of required iterations in Algorithm 2.

We have also implemented a mask-based outlier management method that filters out irregular points that potentially lead to unbounded error estimation. Using the CFD simulation data generated by GE as the example again, for nodes with values of $V_{x}=V_{y}=V_{z}=0$, their decompressed values will be tiny when the error bounds used for retrieval are small. These close-to-zero values, however, could yield loose upper bounds for $V_{\textrm{total}}$ (seeing Theorem 2) despite the small real errors. Accoordingly, in this example, we will use a bit-map to record the position of any data point with $0$ total velocity, and only refactor data points whose values are non-zero.

Despite the fact that the proposed QoI error-control theories and pipeline can be generalized to any progressive compressors that meet the Definition 1, the amount of data retrieved may vary due to the different progressive refactoring and error bounding theories adopted by each compressor. In this section, we review and evaluate the pros and cons of three error-controlled progressive compression algorithms.

Error-controlled compression with multiple snapshots: This type of methods leverages existing error-controlled compressors to compress data using a list of error bounds $\{\epsilon_{i}\}$, ranging from small to large. When an error bound $\epsilon^{*}$ is requested during progressive retrieval, one can choose a snapshot with a minimal $i$ such that $\epsilon_{i}<\epsilon^{*}$. Due to the overlapping information among these snapshots, redundancies may be high when multiple precisions are requested during the progressive retrieval.

Error-controlled delta compression [16]: This type of methods also reduces data into multiple snapshots of multi-precision segments but eliminates the redundancy across the snapshots by compressing the residues (errors between the original and decompressed data) instead of the original data. As a result, it can be more efficient than directly compressing data into multiple snapshots with different error bounds but requires retrieving all first $i$ snapshots ($i$ is the minimal integer such that $\epsilon_{i}<\epsilon^{*}$) when the target error bound is $\epsilon^{*}$.

Progressive compression with bitplane: This type of methods does not require to pre-set error bounds. Instead, it encodes the data using bitplanes such that the data can be retrieved and recomposed on demand. Similar to the error-controlled delta compression, it may not be the most efficient when only a single error bound is requested – directly compressing data using the requested error bound usually generates the smallest data footprints in such cases. PMGARD [15] is the leading technique in this kind, but its performance suffers from loose error control and thus may return more precision fragments than needed.

Below, we evaluate the performance of these three kinds of methods when a set of successively lower relative error bounds (i.e., a series of requests $\{\epsilon_{i}\}$ such that $\epsilon_{i+1}<\epsilon_{i}$) are requested. For the first two categories of progressive compressors, we use $SZ3$ as the underlying error-controlled compressors as it provides the tightest $L^{\infty}$ error bound and thus could yield larger compression ratios than compressors with looser error bound [33, 43]. From now on, we refer to the integration of progression with multiple-snapshot, progression with delta compression, and SZ3 as PSZ3 and PSZ3_delta, respectively. For progression with bitplane, we use PMGARD as the underlying refactor. The evaluations have been conducted using GE’s CFD simulation data and their six QoIs (see Table III for details on datasets).

We evaluated $4$ data fields in Fig. 2. The trend of VelocityY is similar to those of VelocityX and VelocityZ, and thus omitted. We set $\epsilon_{i}=10^{-i}$ for $i=1,2,\cdots 10$ for both PSZ3 and PSZ3_delta to create multiple snapshots (because this setting has a reasonable trade-off between additional storage cost and retrieval efficiency), and requested errors bounds on primary data as $\{\epsilon^{\prime}_{i}\}=0.1*2^{-i}$ for $i=1,2,\cdots 20$. The rate-distortion curves can be interpreted as follows: for any data point $(x,y)$, its value represents the bitrate $x$ (analogous to the percentage of data retrieved) under the requested tolerance $y$. As such, the closer a curve is to the left bottom (left indicates lower retrieve percentage and bottom indicates low error), the better. Notably, PSZ3 has large bitrates under progressive requests, which is expected due to the redundancy among the multiple snapshots. Since PSZ3 and PSZ3_delta rely on pre-set error bounds to produce multi-precision fragments, their bitrate curves exhibit a stair-case pattern: the amount of retrieved data remain constant across several adjacent error bounds and then present a sudden drop as the retrieval moves to the next snapshot. This leads to suboptimal results when the desired error bound is slightly lower than one of the pre-set error bounds. On the contrary, the trends of bitrate in PMGARD are linear with respect to the requested error bounds. Nevertheless, PMGARD constantly generates larger bitrates than those of PSZ3-delta, as the error bounds implemented with the former are looser than the latter.

For PMGARD, the gap between the requested bounds and real errors is mainly caused by the decomposition algorithm, in particular a $L^{2}$ projection, which maps low-level coefficient nodes to high-level nodal nodes employed for data decorrelation. This decomposition algorithm is derived from MGARD, which is specifically designed to provide optimal error control in $L^{2}$ norm, whereas can cause over-pessimistic estimation on $L^{\infty}$. This is further validated by the experimental results shown in Fig 3, which examines the difference among the requested tolerance, estimated upper bound, and actual error measured after progressive retrieval. It shows that while the estimated errors are close to the requested tolerance, the actual errors are far smaller, causing an over-retrieval problem.

We propose to reduce the gap by omitting the $L^{2}$ projection in PMGARD’s decomposition algorithm. This could yield two benefits. First, without the cross-level intervention, the $L^{\infty}$ norm can be accurately estimated through a summation of the maximal error bounds across all levels. Second, because $L^{2}$ projection is time-consuming, removing it can accelerate the refactoring and reconstruction process. We call the new progressive algorithm without $L^{2}$ projection PMGARD-HB as it replaces the orthogonal basis in MGARD with a hierarchical basis from a mathematical point of perspective. As shown in Fig. 3, PMGARD-HB yields much more accurate error estimation than PMGARD, and this leads to constantly lower bitrates in all requested error bounds. We further compare PMGARD-HB with PSZ3 and PSZ3_delta in Fig. 2, where one can observe the advantages in most bitrates.

We conduct all our experiments using the Morgan Compute Cluster (MCC) [44] located at the University of Kentucky 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 memory. Our benchmark datasets are from multiple computational science domains including CFD, cosmology [45], climate [46], and combustion [47], and their specifications are listed in Table III. Note that we use $\{\_\_\}$ in the dimensionality of GE data because its second dimension may have variable sizes. We use the first four datasets for sequential evaluation and the last dataset for measuring data transfer performance in a distributed memory environment. For QoIs, we use Equation (1) – (6) for the two GE datasets, and we test total velocity (i.e., $V_{total}$ in Equation (1), also referred to as VTOT in the rest of the evaluation) for NYX and Hurricane data to demonstrate the generalibility. For S3D data, the data represent the molar concentration of $8$ species associated with $21$ reactions, and their multiplications generate the intermediate variables to derive the rate of progress. In particular, we present $4$ multiplications involved in $2$ reactions. For instance, $x_{0},x_{1},x_{3},x_{4},x_{5}$ used in our evaluation represent species $H_{2},O_{2},H,O,OH$, respectively, and $x_{1}x_{3}$ computes the molar concentration of $O_{2}$ and $H$ in the reaction $H+O_{2}\xrightleftharpoons{}O+OH$. While we convert single-precision NYX and Hurricane data to double-precision for evaluation with small error bounds, our method directly applies to single-precision floating-point data.

We first show that our theory provides guaranteed error control on the derivable QoIs in the evaluated applications. We use PMGARD-HB for demonstration purposes, and the same functionality can be provided by PSZ3 and PSZ3-delta as well. In particular, we present the max estimated QoI errors and actual QoI errors of the proposed method under a progressive set of requested QoI errors for GE data in Figs. 4. It is observed that the actual QoI errors in our method are always smaller than the estimated QoI errors, which are usually close but strictly smaller than the requested QoI errors. This validates our theory, which always provides an upper bound for the QoI error estimations. Different trends are observed for different variables as well. For instance, one can see a gap between the max estimated errors and actual errors in total velocity when the bitrate is low. This is because some decompressed data become close to $0$ when the error bound is high, in which case the estimation of $\sqrt{x}$ generally leads to a loose bound. This situation becomes better when the error bound decreases to a certain threshold, which implies the diminishment of near-zero decompressed values. Furthermore, one can notice a larger gap between the max estimated errors and actual errors in PT than that in the other QoIs, which is reasonable because the estimation in PT is the most complex and involves more relaxation. In addition, the trends in T and C are similar, which is also as expected because their formulas are very similar.

We also present the results of total velocity on NYX and Hurricane data in Fig. 5, as well as four examples for molar concentration multiplications on S3D data in Fig. 6. Similar trends of QoI errors in total velocity are observed in the other two datasets, which demonstrate the generality of our algorithm. Also, our QoI error estimator demonstrates high accuracy on S3D QoIs. This is because these QoIs only involve the multiplications of two variables, which have predictable errors in almost all cases.

We then compare the efficiency of the three progressive approaches in terms of their retrieved data sizes. In particular, we compare their bitrates under one requested QoI error to demonstrate generic cases; requesting progressive error bounds leads to similar results for PSZ3-delta and PMGARD-HB but may negatively impact PSZ3. Similar to the setting in Section V, we choose $\epsilon_{i}=10^{-i}$ for $i=1,2,\cdots 18)$ (used 18 because some datasets such as S3D requires high precision) as the pre-set error bounds for PSZ3 and PSZ3-delta.

We showcase how the proposed method can potentially improve the performance of data retrieval from remote sites using the GE-large data with VTOT as the target QoI. The refactored data is stored at MCC, and the data retrieval request is initiated from the Anvil supercomputer [48] located at Purdue University. The experiment is performed using $96$ cores, each of which will deal with one data block in the GE-large data independently, and the data transfer is performed using the renowned data service software Globus [49]. The data refactoring time for PMGARD-HB, PSZ3, and PSZ3-delta is 2.17 seconds, 5.18 seconds, and 4.67 seconds, respectively, and the total data transfer time is depicted for remote retrieval in Fig. 9. As a baseline, the transfer time of the original data ($3$ variables, 4.67 GB in total) is roughly $11.7$ seconds, as indicated by the dashed line. For progressive approaches, the data transfer time includes the retrieval time, which determines the proper amount of data, and the transfer time, which is the actual time for transmitting them. It is observed that all the progressive approaches lead to less total data transfer time when certain QoI errors can be tolerated. PMGARD-HB and PSZ3-delta exhibit similar performance as they have similar sizes of reduced data in this case, but PMGARD-HB features a shorter data refactoring time as noted above. When compared with the vanilla data transfer with the original data, PMGARD-HB yields $2.02\times$ data transfer performance if the requested QoI error tolerance is 1E-5, because the size of the transferred data is less than 27% of the original one.