Numerical validation of compensated algorithms with stochastic arithmetic

This chapter contains sections titled: Introduction Round-Off Error Propagation Induced by FP Computation The CESTAC Method, a Stochastic Approach for Analyzing Round-Off Error Propagation Stochastic Arithmetic Validation and Implementation of the CESTAC Method Discrete Stochastic Arithmetic (DSA) Taking into Account Data Errors The CADNA Library [36] The Use of the CADNA Library Can the CADNA Library Fail? Conclusion References

Floating-point arithmetic precision is limited in length the IEEE single (respectively double) precision format is 32-bit (respectively 64-bit) long. Extended precision formats can be up to 128-bit long. However some problems require a longer floating-point format, because of round-off errors. Such problems are usually solved in arbitrary precision, but round-off errors still occur and must be controlled. Interval arithmetic has been implemented in arbitrary precision, for instance in the MPFI library. Interval arithmetic provides guaranteed results, but it is not well suited for the validation of huge applications. The CADNA library estimates round-off error propagation using stochastic arithmetic. CADNA has enabled the numerical validation of real-life applications, but it can be used in single precision or in double precision only. In this paper, we present a library called SAM (Stochastic Arithmetic in Multiprecision). It is a multiprecision extension of the classic CADNA library. In SAM (as in CADNA), the arithmetic and relational operators are overloaded in order to be able to deal with stochastic numbers. As a consequence, the use of SAM in a scientific code needs only few modifications. This new library SAM makes it possible to dynamically control the numerical methods used and more particularly to determine the optimal number of iterations in an iterative process. We present some applications of SAM in the numerical validation of chaotic systems modeled by the logistic map.

Numerical Validation of Compensated Summation Algorithms with Stochastic Arithmetic

We present a unified framework for most of the known and a few new evaluation algorithms for multivariate polynomials expressed in a wide variety of bases including the Bernstein-Bézier, multinomial (or Taylor), Lagrange and Newton bases. This unification is achieved by considering evaluation algorithms for multivariate polynomials expressed in terms of L-bases, a class of bases that include the Bernstein-Bézier, multinomial, and a rich subclass of Lagrange and Newton bases. All of the known evaluation algorithms can be generated either by considering up recursive evaluation algorithms for L-bases or by examining change of basis algorithms for L-bases. For polynomials of degree n n in s s variables, the class of up recursive evaluation algorithms includes a parallel up recurrence algorithm with computational complexity O ( n s + 1 ) O(n^{s+1}) , a nested multiplication algorithm with computational complexity O ( n s log ⁡ n ) O(n^s \log n) and a ladder recurrence algorithm with computational complexity O ( n s ) O(n^s) . These algorithms also generate a new generalization of the Aitken-Neville algorithm for evaluation of multivariate polynomials expressed in terms of Lagrange L-bases. The second class of algorithms, based on certain change of basis algorithms between L-bases, include a nested multiplication algorithm with computational complexity O ( n s ) O(n^s) , a divided difference algorithm, a forward difference algorithm, and a Lagrange evaluation algorithm with computational complexities O ( n s ) O(n^s) , O ( n s ) O(n^s) and O ( n ) O(n) per point respectively for the evaluation of multivariate polynomials at several points.

In this chapter, we consider the evaluation of a polynomial function of a single variable. We usually compute the value of an arithmetic function by replacing each arithmetic operation by its corresponding floating-point machine operation (see Section 3.5). Roundoff errors and cancellations sometimes cause the calculated result to be drastically wrong. For similar reasons, a naive interval evaluation of a polynomial may lead to intervals so large as to be practically useless. Roundoff and cancellation errors are especially dangerous if we are evaluating a function close to a root, as we will see in Chapter 9 when we compute verified enclosures of zeros of polynomials.

More than 45 years ago, Dekker proved that it is possible to evaluate the exact error of a floating-point sum with only two additional floating-point operations, provided certain conditions are met. Today the respective algorithm for transforming a sum into its floating-point approximation and the corresponding error is widely referred to as {{,mathrm{FastTwoSum},}}. Besides some assumptions on the floating-point system itself—all of which are satisfied by any binary IEEE 754 standard conform arithmetic, the main practical limitation of {{,mathrm{FastTwoSum},}} is that the summands have to be ordered according to their exponents. In most preceding applications of {{,mathrm{FastTwoSum},}}, however, a more stringent condition is used, namely that the summands have to be sorted according to their absolute value. In remembrance of Dekker’s work, this note reminds the original assumptions for an error-free transformation via{{,mathrm{FastTwoSum},}}. Moreover, we generalize the conditions for arbitrary bases and discuss a possible modification of the {{,mathrm{FastTwoSum},}} algorithm to extend its applicability even further. Subsequently, a range of programs exploiting the wider applicability is presented. This comprises the OnlineExactSum algorithm by Zhu and Hayes, an error-free transformation from a product of three floating-point numbers to a sum of the same number of addends, and an algorithm for accurate summation proposed by Demmel and Hida.

Standard Wilkinson-type error estimates of floating-point algorithms that are solely based on the first or second standard model typically involve a factor γ k := k u / ( 1 − k u ) \gamma _k:=k\mathbf {u}/(1-k\mathbf {u}) , where u \mathbf {u} denotes the relative rounding error unit of a floating-point number system. Using specific properties of floating-point grids it was shown that often γ k \gamma _k can be replaced by k u k\mathbf {u} , and the restriction on k k can be removed. That is true for standard algorithms such as summation, dot product, matrix multiplication, and LU- or Cholesky decomposition. Recently it was shown that, at least for summation and dot product, such results derive without any reference to a floating-point grid. In the current paper we further sharpen the error estimate for summation into k u / ( 1 + k u ) k\mathbf {u}/(1+k\mathbf {u}) , again without any reference to a floating-point grid. Furthermore, an estimate of type h u h \mathbf {u} is shown for sums and dot products that are evaluated using a binary tree of height h h . Both estimates require a mandatory restriction of size 1 / u 1/\mathbf {u} on the number of summands and the height, respectively. Finally, a different kind of error estimate is shown for recursive summation. The discussed bound is sharp, holds true for any number of summands, and is uniformly bounded by 1 1 . The novelty of our approach is twofold. First, rather than using a rounding function, the discussed estimates are based on almost arbitrary perturbations of real operations without any reference to a floating-point grid. As a consequence, the corresponding floating-point error bounds in some base β \beta for rounding to nearest, and partly also for directed roundings, follow as corollaries. Second, in regard to our weak assumptions, the new estimates are sharp. Our main result is sharp for actual realizations of grids floating-point arithmetics are based on. To be more specific, for any feasible problem size, for IEEE 754 754 binary 32 32 as well as binary 64 64 format, there are examples satisfying the given bound with equality.

Systolic algorithm for polynomial interpolation and related problems

The presence of rounding errors is frequently inevitable when performing arithmetic operations in computers due to the use of floating-point number system whose elements have finite precision. Consequently, error-free transformation algorithms are proposed as solutions, and a particular instance of these is the error-free transformation algorithm for real matrix multiplication. Based on the previously mentioned algorithm, this study proposes two error-error free transformation algorithms for complex matrix multiplication. By conducting the numerical experiment, one finds that the combination of both proposed algorithms with an accurate summation algorithm generates the most accurate results, but it also suffers from a notable increase in the computing time ratio. On the other hand, the second proposed algorithm, which is simply another form of the first one with the change of order in summation, combined with the pure floating-point addition produces results that are nearly as accurate as those generated by the former combination. Since these results are achieved with significantly smaller computing time ratio, then the study concludes that the last combination is the best method for complex matrix multiplication in terms of accuracy and computing time efficiency.

Lin, Chung and Han introduced in 2014 the LCH-basis in order to derive FFT-based multipoint evaluation and interpolation algorithms with respect to this polynomial basis. Considering an affine space of n = 2j points, their algorithms require O(n · log2 n) operations in F2r. The LCH-basis has then been extended over finite fields of characteristic p by Lin et al. in 2016 and an n-point evaluation algorithm has been derived for n = pj with complexity O(n · logp n · p). However, the problem of interpolating polynomials represented in such a basis over FPr has not been addressed.

A straightforward discretisation of high-dimensional problems often leads to a curse of dimensions and thus the use of sparsity has become a popular tool. Efficient algorithms like the fast Fourier transform (FFT) have to be customised to these thinner discretisations and we focus on two major topics regarding the Fourier analysis of high-dimensional functions: We present stable and effective algorithms for the fast evaluation and reconstruction of multivariate trigonometric polynomials with frequencies supported on an index set \(\mathcal{I}\subset \mathbb{Z}^{d}\).KeywordsFast Fourier TransformFourier CoefficientTrigonometric PolynomialSparse GridMatrix Vector MultiplicationThese keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

Potentially, data compression techniques may have a broad impact in computing not only by decreasing storage and communication costs, but also by speeding up computation. For many image processing applications, the use of data compression is so pervasive that we can assume the inputs and outputs are in a compressed domain, and it is intriguing to consider doing computations on the data entirely in the compressed domain. We speed up processing by doing computations, including dot product and convolution on vectors and arrays, in a compressed transform domain. To do this, we make use of sophisticated algebraic techniques for evaluation and interpolation of sparse polynomials. We illustrate the basic methodology by applying these techniques to image processing problems, and in particular to speed up the well known splatting algorithm for volume rendering. The splatting algorithm is one of the most efficient of existing high quality volume rendering algorithms; it takes as input three dimensional volume sample data of size N/sup 3/ and outputs an N/spl times/N image in O(N/sup 3/f) time, where f is a parameter known as footprint size (which often is hundreds of pixels in practice). Assuming that the original sample data and the resulting image are stored in the transform domain and can be lossily compressed by a factor /spl rho/ with small error, we show that the rendering of the image can be done entirely in the compressed transform domain in decreased time O(/spl rho/N/sup 3/ log N). Hence we obtain a significant speedup over the splatting algorithm when f/spl Gt//spl rho/ log N.

Development of a Stokes flow solver robust to large viscosity jumps using a Schur complement approach with mixed precision arithmetic

Optimal algorithms for parallel polynomial evaluation

Classical probabilistic rounding error analysis is particularly well suited to stochastic rounding (SR), and it yields strong results when dealing with floating-point algorithms that rely heavily on summation. For many numerical linear algebra algorithms, one can prove probabilistic error bounds that grow as O(nu), where n is the problem size and u is the unit roundoff. These probabilistic bounds are asymptotically tighter than the worst-case ones, which grow as O(nu). For certain classes of algorithms, SR has been shown to be unbiased. However, all these results were derived under the assumption that SR is implemented exactly, which typically requires a number of random bits that is too large to be suitable for practical implementations. We investigate the effect of the number of random bits on the probabilistic rounding error analysis of SR. To this end, we introduce a new rounding mode, limited-precision SR. By taking into account the number r of random bits used, this new rounding mode matches hardware implementations accurately, unlike the ideal SR operator generally used in the literature. We show that this new rounding mode is biased and that the bias is a function of r. As r approaches infinity, however, the bias disappears, and limited-precision SR converges to the ideal, unbiased SR operator. We develop a novel model for probabilistic error analysis of algorithms employing SR. Several numerical examples corroborate our theoretical findings.