{'fraction_non_alphanumeric': 0.0337761925991975, 'fraction_numerical': 0.010753455193936692, 'mean_word_length': 4.543845773603559, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
QDOT: Quantized Dot Product Kernel for Approximate High-Performance Computing

Approximate computing techniques have been successful in reducing computation and power costs in several domains. However, error sensitive applications in high-performance computing are unable to benefit from existing approximate computing strategies that are not developed with guaranteed error bounds. While approximate computing techniques can be developed for individual high-performance computing applications by domain specialists, this often requires additional theoretical analysis and potentially extensive software modification. Hence, the development of low-level error-bounded approximate computing strategies that can be introduced into any high-performance computing application without requiring additional analysis or significant software alterations is desirable. In this paper, we provide a contribution in this direction by proposing a general framework for designing error-bounded approximate computing strategies and apply it to the dot product kernel to develop qdot -- an error-bounded approximate dot product kernel. Following the introduction of qdot, we perform a theoretical analysis that yields a deterministic bound on the relative approximation error introduced by qdot. Empirical tests are performed to illustrate the tightness of the derived error bound and to demonstrate the effectiveness of qdot on a synthetic dataset, as well as two scientific benchmarks -- Conjugate Gradient (CG) and the Power method. In particular, using qdot for the dot products in CG can result in a majority of components being perforated or quantized to half precision without increasing the iteration count required for convergence to the same solution as CG using a double precision dot product.

1. Introduction. Scientific application development is rife with approximation techniques. From spatial and temporal discretization schemes to numerical solution strategies, approximation techniques are quite ubiquitous in scientific computing. In the context of high performance computing (HPC), approximate computing strategies have been typically employed to (a) reduce computational complexity -for example polynomial approximations of functional forms [19] -and (b) make better use of computing resources -for example iterative refinement on Tensor Core Units [6]. For all of these strategies in HPC, the overarching goal is to safeguard the accuracy of the desired quantity of interest, although in some cases, the additional benefit of preserving the convergence properties of the underlying algorithm may be realized [2].
Outside these traditional approximation strategies, other methods have been developed to reduce computational complexity by economizing computation, that is, computing only when necessary. Strategies such as loop perforation [13] are used to skip specific iterations of a loop in a computational kernel in order to reduce the cost. The choice of which iterations to skip may be random or prescribed by some predetermined stride. Another common strategy known as function memoization [11] employs a look-up table to store computed entries of computationally expensive kernels. Subsequent calls to the kernel may then be substituted by retrieving an appropriate entry in the look-up table. While these approximate strategies are more efficient than their exact counterparts, they are not designed with closed form bounds for error introduced by the approximation. This lack of error control precludes them from being considered as is in HPC applications. Significance-aware strategies [20] attempt to mitigate this by accepting an approximation within some prescribed tolerance.
Another class of approximate computing strategies that has gained a lot of interest in recent years is precision-based approximate computing [9,10,17]. Here, the notion of economizing computation is based on computing only on important information, that is, making use of only the important bits. The approach typically follows identifying regions of code or portions of an algorithm where low precision arithmetic or storage may be utilized, subject to some error tolerance. These strategies can lead to efficient mixed-precision algorithms that can efficiently utilize current and emerging HPC hardware.
Existing approximate computing strategies, including the tolerance-based strategies, lack a priori error analysis on their impact on the convergence and the quality of the quantity of interest in the kernels or algorithms they approximate. For most of these methods, this means they need to be deployed as a preprocessing tool to first identify the parameters whose approximation satisfy the required error tolerance, before running the algorithm with these parameters. In this paper, we present a general framework for developing error-bounded approximate computing techniques. We define key terminology that may be used to evaluate other approximate computing strategies, in terms of their effectiveness as an approximation technique and their computational efficiency. We demonstrate this framework by designing an error-bounded quantized dot product algorithm (qdot) as an approximation technique for a dot product kernel. This new approximation strategy can thus be deployed at run time since its error bound is known and can be manipulated at run time. It is worth noting that hardware-based approximate computing strategies have also been developed for specialized hardware such as FPGA's to exploit efficiency by computing in the presence of hardware faults [13,18]. These techniques are beyond the scope of this paper and, as such, are not addressed here.
The paper is organized as follows. In Section 2, we discuss a general framework for designing error-bounded approximate computing kernels and formalize terminologies for measuring the quality of approximate computing kernels. In Section 3, we apply this framework to the dot product kernel to develop qdot. In particular, we provide motivation and pseudocode for qdot together with theoretical analysis capturing the relative approximation error introduced by qdot. In Section 4, we test our implementation of qdot to empirically verify the error bound as well as demonstrate scientific computing applications with dot products that are amenable to extreme levels of quantization and perforation. Specifically, we integrate qdot with CG and Power method implementations to demonstrate the limits of quantization and perforation that can be achieved without affecting the convergence of these applications. We now provide notation used throughout the paper.
1.1. Notation. For an integer n ≥ 1, we write [n] to denote the set {i ∈ N : 1 ≤ i ≤ n}. We denote scalar values by italicized letters, x, and vectors by boldface letters, such as x. Similarly, we denote scalar valued functions by italicized letters, f , and vector-valued functions by boldfaced letters, f . We write to denote the Hadamard, or component-wise, product of two vectors so that (x y) i = x i y i . We write ⊗ to denote floating point multiplication. If S is a subset of indices of x, then we write x S to denote the subvector of x with indices in S. We use the function flexp : R → Z defined by When x is representable as a floating-point number, flexp(x) is the unbiased exponent of x. For x ∈ R n we naturally extend this notation by writing flexp(x) to denote the vector with ith component flexp(x i ). For x ∈ R n write e max (x) and e min (x) to denote the maximum and minimum value of flexp over all components of x. These can be defined mathematically as 2. Error Bounded Approximate Kernels. Suppose we are given a kernel described as a function k : D → R, where the domain and range of k are normed vector spaces. Our goal is to develop a parametric approximation of k for which suitable parameter choices result in bounded error approximations to the output of k. Depending on the complexity of the original kernel k, the level of parameterization required to achieve bounded error in the approximate kernel could easily become nontrivial. Hence, in order to promote usability of any approximate kernel with guaranteed error bounds, This manuscript is for review purposes only.
we also propose that they should be accompanied by a function for identifying parameter choices that ensure the error will be bounded at a specified tolerance.
Let P denote the set of all possible parameter configurations for the approximate kernel, then we can define an approximate kernelk : D × P → R and a parameter selection function Θ : D × R → P such that for any input X ∈ D and approximation error tolerance ε > 0 we have a bound on the relative approximation error given by It is important to formalize the roles thatk and Θ play in achieving approximation in practice. To this end, we establish the following definitions that emphasize the design goals for error-bounded approximate computing kernels and simplify our discussion of them.
Definition 2.1. Let k : D → R be a kernel,k : D × P → R be a parameterized approximation of k, and Θ : D × R → P a parameter selection function fork. Let ε > 0 and X ∈ D be given.
The approximation potential ofk is how close the relative approximation error can be made to the tolerance ε for any choice of parameter θ ∈ P. When k(X) = 0 the approximation potential can be quantified by where θ * is a solution to the optimization problem The approximation effectiveness ofk and Θ is how close the relative approximation error is to the tolerance ε when using the parameter Θ(X, ε) given by the parameter selection function. When k(X) = 0 the approximation effectiveness can be quantified by If the relative approximation error exceeds the tolerance ε with the selected parameter Θ(X, ε), we say that the approximation is ineffective.
The approximation efficiency of Θ is the ratio of the time spent computing the original kernel, k(X), to the time spent selecting parameters Θ(X, ε) and running the original kernel. That is, The approximation efficiency takes on values in the interval (0, 1). Values close to 0 indicate that more time will be required to identify the parameters (inefficient) while values close to 1 indicate that the time spent identifying parameters is fast compared to evaluating the original kernel (efficient).
The approximation speedup ofk is the ratio of the time spent computing the original kernel, k(x), to the time spent computing the approximate kernel. That is,

(2.6)
Note that for this definition we do not include the runtime of Θ(X, ε) as this is to serve as a measure of the maximum possible speedup achievable by the approximate kernelk. This manuscript is for review purposes only.
The design of the parametric approximate kernelk, in particular the complexity of the parameterization, is solely responsible for determining the approximation potential and approximation speedup. The potential is more of a theoretical value, since an optimal parameter, θ * , will often be difficult to obtain for mostk. We can think of the approximation potential to be the best case approximation achievable byk, provided we have unlimited resources to identify the best choice of parameters. The approximation speedup, on the other hand, indicates how fast the computation ofk is compared to the original kernel, k. It does not include the cost of selecting parameters, and one can consider this as the speedup in executingk, if the parameter selection function, Θ, is only used as a preprocessing tool. Values greater than 1 indicate thatk provides a practical speedup over the original kernel.
On the other hand, the design of the parameter selection function Θ determines the approximation effectiveness and approximation efficiency. We can interpret the approximation effectiveness to be the level of approximation actually realized by replacing k(X) withk(X; Θ(X, ε)) in an application. Hence, Θ has the important role of identifying a parameter choice for computing accurately and tightly bounded error in an application containing the kernel k. An approximate kernel will have low approximation efficiency when the asymptotic time to select parameters is much longer than the time to run the original kernel. In order to have a error-bounded kernel with high approximation effectiveness, the user will likely have to sacrifice some level of efficiency. The approximation efficiency defined in (2.5) provides an indication of the overhead in the parameter selection phase required for the design of the approximate kernelk, with respect to the cost of executing the original kernel k. This overhead can be defined as Asymptotic Runtime of Θ(X, ε) Asymptotic Runtime of k(X) so that (2.5) is equivalent to 1 + Asymptotic Runtime of Θ(X,ε) which is bounded in the interval (0, 1). This formulation of approximation efficiency is motivated by the definition of parallel efficiency [4] of a parallel algorithm, which measures the fraction of time for which a processor is doing useful work. Approximation efficiency values close to 1 indicate that the approximate kernel will spend most of the time doing useful work. Efficiency values less than 1/2 indicate that Θ will not be performant enough for runtime-based error-bounded approximation.
The design of the approximate kernel and the parameter selection function go hand-in-hand to form an application that will be effective, efficient, and provide speedup in practice. A well-designed approximate kernel with a poor parameter selection function will result in poor approximation in practice, since the effectiveness of the approximation may not reach the potential. On the other hand, a poorly designed approximate kernel will typically result in ineffective approximation regardless of the design of Θ, since the approximation potential may be limited. Further, it is important to note that an over-parameterized approximate kernel may result in difficulty defining efficient parameter selection functions.
Based on our theoretical formulation of the problem, we summarize our framework for designing error bounded approximate kernels in two phases: 1. Parameter Selection Phase: Define a function or algorithm, Θ, that takes a kernel input X and an error tolerance ε and returns approximate kernel parameters that ensure that the relative error of the approximationk evaluated on the input X is below the tolerance ε. 2. Computation Phase: Use parameters returned by Θ in approximate kernel with input X and error tolerance ε.
It is clear from this general framework that the cost of ensuring bounded error in the approximate kernel is based on the approximation effectiveness and efficiency of Θ.
Before applying this general framework to design an approximate kernel for the dot product, we briefly consider two examples of existing approximate computing techniques and demonstrate how they closely match this framework. First, loop perforation [13] is an approximate computing technique where iterates in a loop are skipped during execution of the code. For a loop of length N, the parameters for this approximate kernel are subsets of [N ] that indicate which iterates of the loop should be skipped. Popular choices for these parameters are selecting random subsets of [N ] or the last k iterates of the loop This manuscript is for review purposes only.
independent of the input X. These choices are often made for efficiency as, in this case, the parameter selection function has asymptotic runtime of O(1), while the original kernel has runtime O(N ). Hence, the parameter selection function is very efficient. However, for many inputs the approximation error will be unbounded, making the kernel ineffective. For applications that are error-sensitive, which is the target of this paper, these types of parameter selection functions are infeasible. Second, approximate computing strategies such as those based on precision tuning and significance-aware methods can be effective, since they are designed to satisfy some error threshold. However, the parameter selection phase of these strategies can be quite expensive, rendering them inefficient for approximate computing at run time. Hence, they are typically run as an offline analysis tool on a representative set of inputs prior to the actual application run. Examples of techniques that make use of such an approach are Precimonious [17], FPTuner [1], and ADAPT [10].
By framing the error-bounded approximate kernels in this two phase approach, it also provides a clear technique for making use of existing unbounded approximate computing methods. In particular, to make unbounded methods feasible in error sensitive applications, we simply need to replace the existing parameter selection function with one that ensures the approximation error is bounded. We now use this general framework to design and test an error bounded approximate kernel for the dot product.
3. QDOT: Quantized Dot Product Kernel. The design of the approximate dot product kernel, qdot, is motivated by our interest in understanding how much approximation can be tolerated by scientific computing applications, without compromising the fidelity of the quantities of interest computed by the application. Enabling easy integration of qdot into existing applications and providing tuneable error bounds, allows for easy evaluation of the approximate kernel in various scientific computing applications. To this end, our primary goal in the design of qdot is to develop an approximate kernel that is effective. If the application can handle large amounts of approximation, then qdot can be used as starting point to develop a more efficient approach.
qdot can be described as a combination of variable precision computing and loop perforation, applied to a full precision dot product kernel. At a high-level, qdot takes the vectors x, y ∈ R n and partitions x and y into sub-vectors of different precisions. Then, for each sub-vector we compute the dot product in its corresponding precision. The dot products from each sub-vector are then accumulated in the same precision as the input, resulting in an approximation of the full-precision dot product kernel. Hence, our approximate algorithm requires two sets of parameters: 1. Disjoint bins {B k } p k=1 specifying indices for p sub-vectors of x and y, where 2. Bin precisions {π k } p k=1 , where π k specifies the precision for the sub-vectors x B k and y B k . For our implementation, the precisions π k can take on IEEE Floating Point Precisions half, single, double, quad as well as perforate. Note that perforate denotes all indices that are to be excluded from the computation of the quantized dot product. In essence, we can consider perforate to be components that have been quantized to a 0-bit representation. Hence, a parameter selection algorithm for qdot should identify bins and the associated precisions for each bin, so that the approximate dot product has bounded error.

Computation Phase.
Given vectors x, y ∈ R n , the objective of qdot is to approximate the value of x y. The design of qdot was motivated by loop perforation but refined to take advantage of other levels of approximation that loop perforation does not account for. For example, given x and y in double precision, loop perforation can only select indices to exclude in the approximation of x y. The use of only perforation as an option to introduce approximation will likely result in a kernel that has limited approximation potential in instances where a small approximation error tolerance ε is required by the application. However, we can improve the approximation potential in these instances by including intermediate levels of approximation between double and perforate by determining component-wise Algorithm 3.1 qdot Computation Phase 1: Inputs: Arrays x, y ∈ R n , parameters θ = {B k , π k } p k=1 , where B k contains indices for sub-vector k and π k denotes the precision of bin B k . 2: Output: Approximate dot product z 3: for k = 1 to p do 4: Initialize sub-vectors with precision π k 5: Compute dot product of sub-vector k at precision π k 6: end for where addition is at precision of x and y Accumulate across sub-vectors products that can be computed at single and half precisions, together with component-wise products that can be perforated. To see this, consider the process of computing all of the component-wise products in the dot product before accumulating the values. Let z = x y, where each component of z is computed at the input precision. Then all components have the same number of mantissa bits, say µ, but potentially varying exponent values given by flexp(z). Now, whenever flexp(z j ) < e max (z), the component z j is a candidate for reduced precision or perforation since the least significant mantissa bits in z j may contribute information outside of the required tolerance. Thus, perforation is suitable only when the dynamic range of the exponents of the entries in z is sufficiently large, or a large tolerance, ε, is acceptable. single and half however, allow the introduction of approximation when the dynamic range is smaller, or a smaller value of ε is required.
The computation phase of qdot computes the dot product of disjoint sub-vectors of the input arrays x and y at some specified level of precision. The bins representing the disjoint sub-vectors, and the precision of each sub-vector, are prescribed by the parameter selection algorithm, and passed as input to the computation phase. Pseudocode for the computation phase of qdot is provided in Algorithm 3.1. In order to simplify the presentation of the error analysis for Algorithm 3.1 and to motivate suitable choices for parameters that will result in bounded error, we provide the following definition: We now establish a bound on the approximation error introduced by Algorithm 3.1.
Theorem 3.2. Suppose x, y are floating point arrays of size n with machine epsilon ε machine . Let {( k , u k ]} k≥1 be a partition of [e min (x y), e max (x y)] and let θ = {B k ,u k (x y), π k } p k=1 denote parameters for Algorithm 3.1. Ifk(x, y; θ) is the output of Algorithm 3.1 then where ε(k) denotes the machine precision for precision π k and M k denotes the cardinality of bin We momentarily consider each inner summation in (3.2). Let ε(k) denote machine epsilon that will be used for the components at the indices in bin B k ,u k (z), for k ∈ [p], where we require each ε(k) ≥ This manuscript is for review purposes only. ε machine . By the fundamental axiom of floating point arithmetic [19] and our choice of ε(k) ≥ ε machine , we have that x j y j − x j ⊗ y j = −x j y j ε, for j ∈ B k ,u k (z) and some |ε| ≤ ε(k). By definition of Algorithm 3.1, it now follows that which is the desired result, (3.1).
Theorem 3.2 establishes that the use of ( , u)-exponent bins for the bin parameters can be used to formulate a closed-form bound on the absolute approximation error. Further, the error is dependent on the size, maximum exponent, and precision of each ( , u)-exponent bin. This allows us to limit the bin parameter search space for the parameter selection algorithm to only ( , u)-exponent bins. By including one additional assumption, we can use this result to establish a relative error bound for Algorithm 3.1.
Ifk(x, y; θ) is the output of Algorithm 3.1 then we have the relative error bound where ε(k) denotes the machine precision for bin precision π k and M k denotes the cardinality of bin Proof. Let E = flexp (x y). By our hypothesis that e max := e max (x y) ≤ E, it follows that 2 −E ≤ 2 −emax which together with the property that |x y| ≥ 2 E yields Combining (3.6) with (3.4) yields the desired result (3.5). Theorem 3.3 provides a way forward for designing the parameter selection phase for qdot. In particular, the parameter selection phase will need to determine the appropriate precisions for each exponent bin that results in the summation in (3.5) being less than the required tolerance. The only practical issue with this result is the required hypothesis e max (x y) ≤ flexp (x y). Of course, this will hold when computing the norm x 2 2 = x x as all of the componentwise products will be nonnegative. For completeness, we provide a relative error bound for when the condition e max (x y) > flexp (x y) is satisfied. However, unsurprisingly, this result requires knowing the value of flexp (x y) to compute the error bound exactly.
is a partition of [e min (x y), e max (x y)] and e max (x y) > E. Then we have the relative error bound where ε(k) denotes the machine precision for precision π k and M k denotes the cardinality of bin This manuscript is for review purposes only.

Parameter Selection Phase.
Since we are interested in approximation effectiveness in this work, the goal of the parameter selection agorithm is to identify the lowest precision that each component of the input arrays can be quantized to, while ensuring the approximation error is below the specified tolerance ε. Based on Algorithm 3.1 and the error analysis in Section 3.1, the key to identifying the appropriate level of precision for each component can be achieved by partitioning the interval [e min (x y), e max (x y)] into p disjoint intervals {( k , u k ]} p k=1 and binning the indices in ( k , u k )-exponent bins. By using the size and maximum exponent of each exponent bin B k ,u k , we can determine the highest level of quantization that can be applied to the indices in each bin without violating the error tolerance requirement.
Suppose that an application requiring the dot product can effectively utilize an approximate dot product algorithm satisfying a relative error tolerance of ε. Then using (3.5) in Theorem 3.3, our goal would be to choose ε(k), for k ∈ [p], such that . This choice of ε(k) corresponds to choosing the precision at which to compute the products for the elements in each bin. Note that for bins with π k = perforate, the value of ε(k) is equal to 1.
If we follow the outline above, we can determine that the precision for each ( , u)-exponent bin is given by the function where the condition involving µ ensures that ε(k) ≥ ε machine is satisfied for each bin. Here, we refer to σ k as the bin score and define it by σ k := log 2 (M k ) + u k − e max − log 2 (ε) + 1 (3.10) for nonempty bins and −∞ for empty bins. It is important to note that the value of µ is dependent on the precision of the input arrays and is defined by and m double , m single , and m half denote the number of mantissa bits in double, single, and half precision, respectively. While these values are known and fixed, we have chosen to use this notation to represent them so that it might be clear how this could generalize to other, perhaps intermediate, levels of precision or alternative number representations [5].
The following theorem establishes the relative error bound achieved by Algorithm 3.1 when (3.9) is used to determine the bin precision parameters {π k } p k=1 . Theorem 3.5. Let x, y be floating point arrays of size n with machine precision ε machine . Suppose that e max := e max (x y) ≤ flexp (x y) and θ = {B ( k ,u k ] (x y), π k } p k=1 , where (3.9) is used with tolerance ε to determine the precisions {π k } p k=1 at which the dot product for each bin is computed. If k(x, y; θ) is the output of Algorithm 3.1 then there exists τ = O(ε) such that x y −k(x, y; θ) ≤ τ |x y| . (3.12) This manuscript is for review purposes only.
Proof. Let b k = log 2 (M k ), for all k ∈ [p]. By our hypotheses that e max ≤ flexp (x y), the result in (3.5) from Theorem 3.3 holds. Since the goal is to bound the relative approximation error of x y by ε, we show that the choices for bin precisions given in (3.9) provides the desired result up to a scalar multiple. Let k ∈ [p], then we require a choice for ε(k) such that Since M k 2 u k −emax+1 ε(k) ≤ 2 b k +u k −emax+1 ε(k) and ε ≥ 2 log 2 (ε) , we can ensure (3.13) holds by enforcing which, when rearranged, yields the bound ε(k) Reordering and labeling components allows for a meaningful interpretation of the bound on µ k as As µ k denotes the number of mantissa bits required to achieve the desired bound, the inequality in (3.15) allows us to determine the minimum level of precision the components in bin B k ,u k (x y) can be computed at while satisfying the desired bound. When the value b k + u k − e max − log 2 (ε) + 1 falls between the number of mantissa bits for two different IEEE Floating Point representations, it is necessary to choose the precision with more mantissa bits to ensure the desired relative error bound is satisfied. If b k + u k − e max − log 2 (ε) + 1 ≤ 0 then the components in bin B k ,u k (x y) are not required to compute the approximate dot product at the desired level of accuracy, which results in assigning a precision of perforate to bin k. Following these observations yields the definition of the Precision function in (3.9) and the bin score value in (3.10). Now using (3.9), we have that M k 2 u k −emax+1 ε(k) ≤ ε, for all k ∈ [p], and hence, where we write N bins to denote the number of nonempty bins and note that N bins ≤ p. Finally, substituting (3.16) into (3.5) yields |x y −k(x, y; θ)| ≤ N bins ε |x y|, the desired result.
From a purely theoretical standpoint, the entire parameter selection algorithm can be summarized by using Definition 3.1 to compute the exponent bins, then using (3.9) to compute the precision for each bin. Pseudocode for the parameter selection algorithm is given in Algorithm 3.2.
Line 3 of Algorithm 3.2 computes the array of exponents flexp(x) + flexp(y), which is the key piece of information required to bound the relative approximation error. In our implementation of the parameter selection phase, this step is accomplished by using a bitmask and bitshift to extract the unbiased floating point exponent of x i and y i , and saving their sum in an array of integers. For inputs with precision up to quad, these sums can be computed using a short int since the quad exponent is 15 bits long and an unsigned short int has a capacity of 16 bits. Furthermore, by working in the exponent space, we can avoid introducing additional floating point operations in the parameter selection algorithm. In fact, nearly all of the operations required by our exponent preprocessing algorithm involve only integers, in particular short int values, and not floating points. This step has a complexity of O(n).
Algorithm 3.2 requires that the exponent array, e, is sorted. In our implementation, we make use of the counting sort algorithm based on two observations. Firstly, the counting sort has a worst-case Algorithm 3.2 qdot Parameter Selection Algorithm 1: Inputs: Arrays x, y ∈ R n , number of mantissa bits in precision of inputs µ, approximation error tolerance ε, and strategy for partitioning exponent range. 2: Outputs: Bins and corresponding Bin Precisions θ ← {B k ,u k , π k } p k=1 3: e ← flexp(x) + flexp(y) Exponent preprocessing 4: e min ← min{e i : 1 ≤ i ≤ n}; e max ← max{e i : 1 ≤ i ≤ n} Identify min and max exponents 5: e ← sort e Sorted bin initialization 6: B k ,u k ← {j ∈ [n] : k < e j ≤ u k } Identify components for each exponent bin 7: σ k ← log 2 (M k ) + u k − e max − log 2 (ε) + 1 Compute bin scores from (3.10) 8: π k ← Precision (B k ,u k (x y)) Determine bin precision using (3.9) runtime of O(n + c), where n is the number of elements being sorted and c is the number of distinct values that the n elements can have. In Algorithm 3.2, the number of distinct values being counted corresponds to the exponent range, which is given by c = e max − e min + 1. Hence, the worst-case cost for using counting sort is O(n + e max − e min + 1). Secondly, the size of the exponent range is bounded fairly tightly regardless of the precision used. For example, input arrays with double precision have exponents in [−1023, 1024] which yields that e max − e min + 1 ≤ 2048 = 2 11 . Further, input arrays with quad precision satisfy e max − e min + 1 ≤ 32768 = 2 15 . Compared to sorting algorithms with runtime O(n log 2 n), counting sort should be theoretically faster for problems of dimension n = O(1e3) and larger, regardless of input precision. Hence, the runtime and memory complexity for this step is O(n + e max − e min + 1).
In line 6 of Algorithm 3.2, we split the original vector into sub-vectors using exponent bins that depend on a partition of the interval [e min , e max ]. To simplify the parameter selection algorithm, instead of explicitly requiring a partition of this range, we specify a strategy that is used to determine partition parameters k and u k during the parameter selection algorithm. After considering various strategies, we follow a straightforward approach that emphasizes approximation effectiveness: • Exact binning: The interval [e min , e max ] is partitioned into singleton intervals defined by where u k = e min + k − 1.
In addition to having the lowest asymptotic runtime and memory complexity compared to other binning strategies we considered, exact binning also helps our implementation avoid overflow and underflow issues when converting input data from higher precisions to lower precisions (e.g., double to half).
To avoid running into these issues, we use the following technique outlined here for bin B k ,u k . First, the data for B k ,u k is scaled so that the unbiased floating point exponent is 0. The scaled data is then copied into an array with precision π k . The dot product is then computed for the components in bin B k ,u k at precision π k . The resulting dot product is then converted back to the original input precision. Finally, by using exact binning, we can undo the scaling by multiplying the resulting dot product for bin B k ,u k by 2 u k . The computational and memory complexity of exact binning is O(n + e max − e min + 1).
The final steps in lines 7 and 8 of Algorithm 3.2 are used to identify the precision of each bin. For p bins, where p ≤ n, this step has a computational and memory cost of O(p). Combining the cost for each step yields the total computational and memory cost of O(n + e max − e min + 1). Since e max − e min + 1 is bounded by 2 11 and 2 15 for double and quad, respectively, it follows that for even moderately sized problems the computational and memory complexity for Algorithm 3.2 are effectively O(n). It is important to emphasize that lines 3 -6 and line 8 make use of short ints. Furthermore, even though line 7 makes use of log 2 , in practice, we found the number of bins, p, to be relatively small (e.g., between 15 and 40).
Before proceeding to the full qdot algorithm, we consider a practical issue that may arise when using this error-bounded approximate kernel. Specifically, it is possible that the parameter selection phase will recommend that all components should remain at the same precision as the inputs. In this This manuscript is for review purposes only.

Algorithm 3.3 qdot -Quantized Dot Product Kernel
1: Inputs: Arrays x, y ∈ R n , and error tolerance ε 2: Output: Approximate dot product z 3: --Parameter Selection Phase --4: Use Algorithm 3.2 with input (x, y), tolerance ε, and exact binning to identify parameters θ 5: --Computation Phase --6: Use Algorithm 3.1 with input (x, y) and parameters θ to compute z case, our approximate kernel will be identical to the original kernel but will incur the overhead cost of the full parameter selection phase. We can establish necessary conditions for when the parameter selection algorithm will recommend using the input precision for all components based only on the values of e min and e max to perform early termination of the quantization.
Theorem 3.6. Let x and y be arrays of floating point values and letμ denote the number of mantissa bits for the level of precision directly below the precision of the inputs x and y and let ε be the tolerance for qdot. If e max − e min ≤ − log 2 (ε) −μ then the level of precision specified for all bins by (3.9) will be the same as the precision of x and y.
Proof. By our hypothesis that e max − e min ≤ − log 2 (ε) −μ we haveμ ≤ − log 2 (ε) − e max + e min . Combining (3.10) together with the fact that log 2 (M k ) ≥ 0 and u k ≥ e min for all k ∈ [p], it now follows thatμ ≤ − log 2 (ε) − e max + e min ≤ − log 2 (ε) − e max + log 2 (M k ) + u k ≤ σ k , for all k ∈ [p]. Sincê µ is the number of mantissa bits for the level of precision directly below the precision of x and y, the Precision function in (3.9) yields that the precision of every bin will be the same as the input.

qdot Pseudocode and Implementation Details.
Putting all the components together, pseudocode for our error-bounded approximate dot product kernel framework, qdot, is presented in Algorithm 3.3. qdot is implemented using C++ and CUDA. The parameter selection and computation phases both run on a GPU to take advantage of parallelism and half precision in hardware. Many steps of the parameter selection phase, as well as the entire computation phase, are amenable to parallelization.
In the implementation of the parameter selection phase, for components x i and y i , the short int exponent values are summed on each thread to determine e i , where each index i ∈ [n] is operated in parallel from the other indices. The values of e min and e max are then computed using a reduction over the array of exponents e. Here, we use thrust::minmax element from the CUDA Thrust C++ template library [16] to perform this reduction. Next, we perform counting sort of e so that the exponent values are in non-decreasing order. We implemented our own version of counting sort in CUDA to take advantage of parallelism. In our implementation of the counting sort algorithm, we save the cumulative counts of each value in the range [e min , e max ], which makes it straightforward to define the bins. All of the remaining steps of the parameter selection algorithm are applied to each bin independently, and hence, are easily parallelizable on the GPU.
While the computation phase is theoretically straightforward, we had to carefully handle computation in lower precisions to avoid overflow and underflow. As noted previously in Section 3.2, this requires scaling the input data prior to performing the quantization, and re-scaling the result back after computing dot product of the individual bins. Each bin dot product is computed using cublas functions from the CUDA Basic Linear Algebra Subroutine library [14]. In particular, as the higher precision values, such as double, do not need to be scaled in our experiments, we compute all of the bins identified as double using a single call to cublasDdot. The lower precision bins, half and single, need to be scaled individually so we use strided and batched cublas functions to efficiently compute the dot products for each bin. This way, they can be re-scaled individually before accumulating. For single bins, we use cublasSgemmStridedBatched, and for half bins, we use cublasGemmStridedBatchedEx. To avoid overflow when computing the dot product for large bins in half precision, we make use of Tensor Core Units [15,22] when calling the function cublasGemmStridedBatchedEx. This is a hardware efficient operation for performing the products in half precision while performing the accumulation in single which avoids overflow. When the Tensor Core Units were not used, we observed that overflow occurred in some cases where the input arrays were large in dimension.

Numerical Experiments.
In this section, we perform numerical experiments using qdot to verify the theoretical error bound and to test the level of approximation that can be used in scientific computing applications. The first set of experiments is performed using synthetic data and serves to empirically verify the error bound and compute the approximation effectiveness, efficiency, and speedup of qdot. Then we integrate qdot into two scientific computing applications to investigate the level of approximation that can be introduced into the dot product kernel without compromising the fidelity of each application. The applications we consider are the Conjugate Gradient (CG) method for solving linear systems, and the Power Method for computing extreme eigenvalues of linear systems.

Tightness of Error Bound.
In the first set of experiments, we used qdot to approximate a double precision dot product for arrays of varying dimension sampled from two different distributions. Since the level of approximation that can be introduced is dependent on the exponent range [e min , e max ], the floating point exponents for the components of the arrays are sampled from uniform or normal distributions, parameterized by a single value t: The theoretical analysis, suggests that larger exponent ranges [e min , e max ] provides the potential for more components to be quantized to lower precisions, or be perforated even when small qdot tolerance values, ε, are used.

Distribution
In the first experiment, we vary the problem dimension while keeping the qdot tolerance fixed at ε = 1e-16, in order to investigate the tightness of the error bound. For each problem dimension, we sample 1900 array pairs from distribution A and 1500 array pairs from distribution B. These pairs are sampled with varying values for the distribution parameter t in order to capture a wide range of input arrays with each distribution. For distribution A, t is taken to be values from the set {5i : 2 ≤ i ≤ 20} with 100 pairs sampled for each t. For distribution B, t is taken to be values from the set {2i : 1 ≤ i ≤ 15} with 100 pairs sampled for each t. For each array pair, we use qdot to compute an approximate dot product value and the cublas function cublasDdot to compute the double precision dot product. The relative approximation error, denoted by relerr in the figures that follow, is equal to the quantity |cublasDdot(x, y) − qdot(x, y; ε)|/|cublasDdot(x, y)|. The relative approximation error bound for qdot follows from (3.12). Since the number of non-empty bins, N bins , is known prior to determining the bin precision, the bin score formula in (3.10) is slightly adjusted to use ε/N bins in place of ε. This is done to ensure that the relative approximation error bound is satisfied for the desired ε. To capture the tightness of the qdot error bound, we compute the difference between the qdot error bound, denoted qdot bound , and the relative approximation error of qdot. For each problem dimension from the set {1e2, 1e3, 1e4, 1e5, 1e6, 1e7}, we generated bar and whisker plots of this difference over all of the samples in Figure 1. We find that data where the exponents are sampled from a normal distribution yields instances where the error is tightly bounded by the error bound even with increasingly large array sizes. On the other hand, the data in which the floating point exponents are uniformly distributed is not captured tightly by the error bound as the problem size increases.
In the following experiments, we evaluate approximation effectiveness, efficiency and speedup on our sampled data problem. We vary the qdot tolerance ε, from 1e-15 to 1e-1 successively by a factor of 10, for data sampled from distributions A and B, each using three different parameter values, t ∈ {5, 9, 13}. The experiments were carried out for problems of dimension {1e5, 1e6, 1e7, 1e8}, and This manuscript is for review purposes only.  For each problem dimension and distribution, we provide a box and whisker plot for the value (qdot bound − relerr) after using qdot to approximate the dot product and cublasDdot to compute the double precision dot product for each sampled array.
for each qdot tolerance value, we sample 50 pairs of arrays. Here, the results were similar for both distributions so we omit results for A. Figure 2 illustrates the results for approximation effectiveness, obtained by computing the relative error, relerr, between qdot and the double precision dot product. For the plot of each problem dimension, the solid line represents the mean and the shaded region of the same color represents the min and max corresponding to that qdot tolerance. In each of these plots, we include a dashed line denoting the qdot tolerance. The approximation effectiveness is then the distance, or gap, between the qdot tolerance and these plots. While the approximation effectiveness appears to be somewhat dependent on the array dimension, particularly for small t, we find that the distribution parameter, t, plays a larger role in the approximation effectiveness. In particular, when sampling from a wide enough exponent range, the approximation effectiveness approaches the optimal value and is essentially the same for all of the array dimensions we tested. In evaluating the approximation efficiency, we observed that for all the test cases considered in these experiments, the computed approximation efficiency was essentially constant at roughly 1 40 = 0.025. As a result, we omit the figures here. This value of approximation efficiency indicates that the parameter selection algorithm is expensive compared to the double precision dot product, and precludes its practical use in applications at runtime. However, there is potential to improve the efficiency of qdot by designing alternative parameter selection algorithms. The most costly step of the current parameter selection algorithm is the sorting step in line 5 of Algorithm 3.2. Designing a parameter selection algorithm that either reduces this cost or does not rely on sorting could result in a more efficient approximation technique, albeit at some cost to its approximation effectiveness.

qdot with Conjugate Gradient.
Conjugate Gradient (CG) [3,8] is a popular iterative technique for solving the linear system Ax = b, given a positive definite matrix A ∈ R n×n and a vector b ∈ R n . Recent works have focused on introducing approximation into linear solvers, such as CG, without compromising the performance or convergence properties of the solver [10]. In the following experiments, we evaluate the impact of introducing approximation into CG by integrating qdot in a double precision CG application. Excluding the matrix vector product, each iteration of CG requires the computation of two dot products. By replacing these dot products with qdot, we can tune the amount of approximation incorporated into the resulting approximate linear solver kernel by varying the value of the qdot error tolerance. Algorithm 4.1 presents pseudocode for our Approximate Conjugate Gradient (ACG) kernel, specifying the exact operations computed using qdot. The implementation of Algorithm 4.1 was obtained by using our implementation of qdot in place of the double precision dot Algorithm 4.1 ACG: Approximate Conjugate Gradient Kernel with qdot 1: Inputs: Array b ∈ R n , symmetric positive definite matrix A ∈ R n×n , initial guess for solution x 0 ∈ R n , CG error tolerance τ , qdot error tolerance ε.
In our experiments with ACG, we solve problems of varying dimension with varying qdot and CG tolerances and visualize how the precision utilized by qdot affects the number of iterations required for convergence. We use the HPCCG application code to generate four different problems: Two 2D problems with dimension 100 × 100 × 1 and 1000 × 1000 × 1 and two 3D problems with dimension 100 × 100 × 10 and 1000 × 1000 × 10. Note that the matrix A generated for each problem is of dimension n × n, where n is equal to the product of the problem dimensions. Hence, solving these problems using CG requires computing dot products for arrays of dimension 1e4 and 1e6 for the 2D problems and 1e5 and 1e7 for the 3D problems. For the examples evaluated here, a preconditioner was not needed as the matrices were relatively well-conditioned. In all of these experiments, we used a fixed linear solver convergence tolerance of τ = 1e-8. These experiments seek to ascertain how much approximation can be introduced into the CG algorithm, using qdot, while preserving its convergence properties. Thus, we vary the qdot tolerance successively by a factor of 10 from 1e-16 to 1e3, to identify the largest tolerance that resulted in ACG converging in the same number of iterations as the double precision CG solver, which we denote simply by CG.
By choosing the largest qdot tolerance, ε, yielding the same number of iterations as CG, the resulting ACG solver utilizes the maximum amount of approximation that qdot can introduce, without impacting the convergence of the original CG solver. For the 2D problems of dimension 100 × 100 × 1 and 1000 × 1000 × 1, this corresponded to ε = 1e0 and ε = 1e3, respectively. For both 3D problems, of dimension 100 × 100 × 10 and 1000 × 1000 × 10, this corresponded ε 1e2. We note that qdot tolerance values exceeding 1.0 are not unusual since the error bound for qdot captures the worst-case scenario error for all arrays of size n. Thus, in some cases, the qdot tolerance could be loose, and yet provide an acceptable amount of approximation.
To illustrate the approximation introduced by qdot, we plot a breakdown of the percentage of each precision used by qdot at each iteration of Algorithm 4.1. Figure 4(a) -(d) presents results for computing r T r and p T Ap for the 2D problems and Figure 4(e) -(h) presents corresponding results for the 3D problems. For the 2D problems, the results indicate that a significant fraction of the dot product computation can be performed in half precision for the smaller dimension problem. For the larger dimension 2D problem, the results suggest a larger fraction could be perforate. For the 3D problems, we observe that a majority of the computation can be performed in half precision or perforate. As the problem size increased, we observe the need to compute more components of the dot product in single precision. We note that while these results are not a rigorous analysis of the Algorithm 4.2 APM: Approximate Power Method Kernel with qdot 1: Inputs: Array b ∈ R n , matrix A ∈ R n×n , initial guess for eigenvector x 0 ∈ R n with x 0 = 1, Power Method error tolerance τ , qdot error tolerance ε. 2: Output: Approximate eigenvalue λ ∈ R and eigenvector x ∈ R n for A. 3: Initialize λ and index 4: while |λ k − λ k−1 | > τ do 5: x k+1 ← Ax k Update guess for eigenvector 6: c ← qdot(x k+1 , x k+1 , ε) Approximate dot product required for normalization 7: x k+1 ← x k+1 / √ c Normalize guess for eigenvector 8: λ k+1 ← qdot(x k , x k+1 , ε)/c Update guess for eigenvalue 9: k ← k + 1 Increment k 10: end while effect of approximating operations in CG, they provide insights into how approximating the dot product operation may impact CG. In particular, they underscore the importance of quantization in preserving the convergence properties of CG. For example, it is easy to see from the results that for larger 3D problems, one cannot guarantee that the convergence rate of CG would be preserved by using solely half precision dot products.
To complete the CG experiments, we provide time lapses of the precision corresponding to each component of the two qdot calls in ACG, on the 2D problem with dimension 100 × 100 × 1. Figure 5a provides results for computing r k r k in line 10 of Algorithm 4.1, and Figure 5b illustrates computing p k q k line 7 of Algorithm 4.1. These figures reveal some structure in the precision distribution of the data used to compute the approximate dot product, as the CG iterations progress. For this particular linear system under consideration, we observe bands of perforate components that increase in size from the boundary of the mesh, as the iteration count increases. As ACG approaches the final iteration some half precision begins to mix with the bands perforated components.

qdot with Power Method.
For the final experiment, we integrate qdot into the Power Method [12,21] to compute the eigenvector-eigenvalue pair corresponding to the largest eigenvalue of a randomly generated graph Laplacian matrix. Like the CG method, the Power Method is an iterative technique where two dot products occur at each iteration. Hence, we experiment with an Approximate Power Method (APM) that is designed by replacing the dot products in the Power Method with qdot. For clarity, we provide pseudocode for APM in Algorithm 4.2.
We fix the convergence tolerance of the power method to τ = 1e-6 and perform a maximum of 300 iterations. Figures 6 and 7 present various results for problems of dimension 1e4 and 1e5, respectively.
For each problem dimension, we vary the qdot ε, and plot the corresponding number of iterations required for convergence. This allows us to identify the highest level of approximation that can be introduced for the same convergence rate as the double precision Power Method application, which we will denote simply by PM. Figures 6a and 7a illustrate the results. In Figures 6b and 7b, we plot the corresponding absolute error between the eigenvalue estimate computed by APM and the eigenvalue estimate computed by PM. This illustrates the effect of ε on the quality of the computed eigenvalue using APM. The results generally indicate that introducing more approximation by increasing the value of ε can lead to differences between the computed eigenvalue using APM and PM. While this behavior generally makes sense, it is also understandable since the stopping criterion used in APM and PM is based on the difference between the eigenvalues computed at successive iterates of the respective algorithms (see Algorithm 4.2). Figures 7c -7d present the distribution of the percentage of components at each precision, corresponding to ε = 1e-7. This is the largest qdot tolerance that resulted in an eigenvalue estimate with an absolute error within the accepted tolerance τ = 1e-6. For these plots, note This manuscript is for review purposes only.   4: qdot Precision at Each Iteration of ACG. We plot the percentage of components at each precision within qdot at each iteration of ACG when using the largest qdot tolerance for which ACG converged in the same number of iterations as CG. Plots in (a) -(d) correspond to 2D problems while (e) -(h) correspond to 3D problems. In any iteration where zero components were assigned a precision, that precision is not plotted at that iteration (to differentiate from very small percentages).  that we only show the nonzero percentages of the precisions that are represented in the components. Interestingly, the results indicate that a very low percentage of components are needed in double precision after the initial iterations.

Conclusion.
In this paper, we formalized a general framework for designing error bounded approximate computing kernels, and applied this framework to the dot product kernel to design qdot. We theoretically proved and empirically demonstrated that qdot bounds the relative error introduced by the approximation. Numerical experiments on the Conjugate Gradient and Power Method algorithms demonstrate that high levels of approximation can be introduced into these algorithms without degrading their performance.
While qdot was demonstrated to have good approximation effectiveness, achieving this effectiveness comes at the cost of a parameter selection phase that limits the approximation efficiency in practice. However, for applications that can tolerate high levels of approximation in the dot product, such as those presented in this paper, one can design a more efficient strategy at some cost to the effectiveness. One strategy is to consider a blocked approach that targets the least efficient step in our implementation -sorting the array of exponent values -to improve efficiency. Future work could explore finding a better balance between approximation effectiveness and efficiency so that qdot could be practically used at runtime. If successful, this could lead to an approximate strategy for the matrix-vector multiplication kernel that balances effectiveness and efficiency. Furthermore, the formalized general framework presented in this paper could be used to design other error-bounded approximate computing strategies for kernels other than the dot product.