Testing the Sharpness of Known Error Bounds on the Fast Fourier Transform

The computation of Fast Fourier Transforms (FFTs) in floating-point arithmetic is inexact due to roundings, and for some applications it can prove very useful to know a tight bound on the final error. Although it can be almost attained by specifically built input values, the best known error bound for the Cooley-Tukey FFT seems to be much larger than most actually obtained errors. Also, interval arithmetic can be used to compute a bound on the error committed with a given set of input values, but it is in general considered hampered with large overestimation. We report results of intensive computations to test the two approaches, in order to estimate the numerical performance of state-of-the-art bounds. Surprisingly enough, we observe that while interval arithmetic-based bounds are overestimated, they remain, in our computations, tighter than general known bounds.


I. INTRODUCTION
The Fast Fourier Transform (FFT) is a ubiquitous tool in many areas of science and engineering.Although it is known [1] to be a fairly accurate algorithm, several applications require a very careful error analysis.A typical example we are interested in is the multiplication of very large integers, or of polynomials of very large degree with integer coefficients, using the Schönhage-Strassen algorithm introduced in [2] or a variant of it: since we know that the exact final values (either digits or polynomial coefficients) are integers, if we are certain that the final error of the (inherently inaccurate) calculation of them in floating-point arithmetic through Fourier transforms is less than 1/2, we can retrieve them.In such cases, the classical error bounds of the form α u + O(u 2 ) that are omnipresent in numerical analysis (see later on for the definition of u) do not suffice: we need to be certain that the error is not larger than a given threshold.For the very same reason, probabilistic error bounds are not considered here, although they are useful in other contexts.
To that purpose, we consider two possible approaches: 1) use global, a priori error bounds, i.e., valid for all possible input values, such as the ones given in [3], and then compute FFTs in conventional floating-point arithmetic (if the global bound is less than the threshold); 2) by computing the FFTs in interval arithmetic (IA) [4], [5], obtain on the fly a local error bound (i.e., valid only for the considered input value).
The first approach is simple and fast.Its main drawback is that by nature a global, uniform, error bound is in general much larger than the actual error we encounter in practice: we may therefore (when the bound is larger than the threshold) erroneously consider, for a given input value, that the floatingpoint calculation will not be accurate enough.The second approach is by nature significantly slower (the ratio depends much on the environment: assuming the precision p is the same-see definition below-interval arithmetic is between around 2.5 and 10 times slower [6], [7] than floating-point arithmetic), and there are two antagonistic effects to consider: • since the computed error bound is local there is hope of obtaining a bound significantly smaller than the global bound; • on the other hand, the tendency of IA to overestimate the size of intervals may result in computed local error bounds significantly larger [8] than the actual error one would obtain by using floating-point arithmetic with the same input values.To choose which of the two approaches should be used for a given problem, we need to know which of these antagonistic effects prevails over the other one, i.e., we need to compare the best known global bounds with local bounds estimated with interval arithmetic.This is the major purpose of this paper.
The literature has a long history [3] of successively-refined global bounds for the Cooley-Tukey FFT algorithm [9] in fixed-point or floating-point arithmetic, that started [10] just one year after Cooley and Tukey's seminal paper.An excellent recent overview of Fourier-related transforms is given in [11].
In [3], the authors present the state-of-the-art global error bounds.They also build a bad case, i.e., an input vector for the FFT for which particularly large errors are obtained, showing that their global bounds cannot be refined much more.
In the following, we compare actually attained errors of the Cooley-Tukey algorithm in binary floating-point arithmetic with global and local bounds obtained through IA, by means of intensive numerical computations.

A. Floating-point arithmetic
We assume that we compute Fourier Transforms in binary, precision-p, floating-point (FP) arithmetic, assuming an unbounded exponent range (hence our results apply to usual FP arithmetics such as the ones specified by the IEEE 754 Standard [12] provided that underflow and overflow do not occur).Throughout this paper, a FP number is a number of the form x = M 2 p−1 • 2 e , where |M | ≤ 2 p − 1 is an integer.Typical values of p are 24 (for binary32/single-precision arithmetic) and 53 (for binary64/double-precision arithmetic).Since the exact sum and product of two FP numbers are not, in general, FP numbers, they must be rounded.Assuming rounded-to-nearest operations (which is the default in IEEE-754 arithmetic), the implementations ⊕ and ⊗ of the addition and multiplication are such that (barring underflow/overflow) for any FP numbers x and y, where u = 2 −p .Eq. ( 1) is sometimes called the "standard model" of FP arithmetic, and a large part of numerical error analysis is based on it [1].In the following, F is the set of the FP numbers, and G is the set of the complex numbers whose real and imaginary parts are FP numbers.

B. Discrete Fourier Transform and Notations
For any positive integer , the discrete Fourier transform takes X = (x 0 , . . ., x N −1 ) ∈ C N and uses the N -th roots of unity An FFT algorithm takes an input vector X ∈ G N and uses FP approximations of ω k N and (rounded) arithmetic operations to return a FP approximation of Y , noted Ŷ ∈ G N .Our main concern is to bound the error | Ŷ − Y |.As there are several possible variants of the FFT algorithm, we consider the one presented in Fig. 1.We assume that N is a power of two.
Function reverse(n,j) for j from 0 to 2 n − 1 do / * Butterfly permutation * / y[j] = x[reverse(n,j)]; for k from 1 to n do y = OneStep(y,k,n); return y Fig. 1.Pseudocode for the radix-2 FFT algorithm, extracted from [3], with permission.Here, reverse(n, j) is the n-bit number whose binary representation is the mirror image of the representation of j, and Nblocs is the number of independent order-2 k FFTs.

C. Different measures of FFT errors
The input and output values of the FFT belong to C N .In that domain, different norms can be used to express distances between vectors (and therefore errors).It is therefore important to precisely establish which norms are used, to avoid confusions.Unfortunately, the "natural norm" for Fourier transforms (i.e., the one for which calculations are easiest) is not the norm that most interests us for our applications.The natural norm for FFT is the 2-norm, used to derive the 2-norm relative error: where However, the 2-norm relative error does not tell us much about the accuracy of a particular coefficient.To get that, we prefer using the following component-wise ∞-norm inputscaled absolute error: where

The choice of scaling by ∥X∥ ⊥
∞ comes from the fact that it makes e ⊥ ∞ comparable to e 2 and, since ∥X∥ ⊥ ∞ is known exactly, an absolute error bound can easily be obtained from (2).In the literature, many results are on e 2 but few on e ⊥ ∞ .Hopefully, one can use: to compare e ⊥ ∞ with e 2 .Note that for large integer multiplication applications, ∥X∥ ⊥ ∞ cannot be arbitrarily large: it is straightforwardly bounded using the value of the largest allowed digit.

III. GLOBAL BOUNDS
In this section, we remind some of the results presented in [3].There are two main sources of errors in an FP implementation of the Cooley-Tukey algorithm: the roundings in the arithmetic operations, and the fact that, except in trivial cases (namely, ±1 and ±i), the roots of unity are not exactly representable in FP arithmetic.

A. Rounding in arithmetic operations
The FFT algorithm uses complex additions and multiplications.From Eq. ( 1), one easily deduces that for x, y ∈ G, the relative error of the floating-point addition ⊕ satisfies Complex multiplication is less simple because there are two implementations to consider: one uses the fused multiply-add (FMA) operation (⊗ FMA ), while the other does not (⊗ MUL ).They are described in [3].For any x, y ∈ G, the errors bounds, proven in [13] and [14] respectively, are In the following, we avoid to specify which implementation is used, and we simply write this bound ρ.

B. Rounded values of roots of unity
FFT algorithms use precision-p FP approximations ωk N of the roots of unity.We note δ n the maximum absolute error of the FP approximations of the order-2 n roots of unity:

C. Global error bound on the Cooley-Tukey algorithm
The state-of-the-art 2-norm global relative error bound is provided by [3].The authors of that paper wrote (with notations adapted): Let Ŷ be the computed 2 n -point FFT of X ∈ C2 n , and let Y be the exact value.Then where g 1 = g 2 = 0, g j = δ j + ρ (1 + δ j ) for j ≥ 3.An ∞-norm input-scaled absolute error can be deduced using Eq. ( 3):

D. A bad case
It is not known if ( 4) is tight, i.e., if there exists a worst case X for which b n nearly equals the actual FP error.However, a bad case, i.e., an input X for which the FP error is particularly large, is built in [3], and the corresponding error is One can see in Fig. 2 that b n and w n have a similar asymptotic behavior with an asymptotic ratio b n /w n of 8.

IV. LOCAL, IA-BASED, BOUNDS
To obtain local bounds, we use a straightforward implementation of the FFT algorithm using interval arithmetic (IA).In IA, a number is represented by an interval z = [z, z].Arithmetic operations in IA, are defined such that the intervals always contain the exact result.In the FFT, we replace the input X by a vector of intervals, the roots ω k N by intervals, and all intermediate operations by their IA equivalents, then we obtain an interval result that contains the exact result: Y ∈ Y.
In computer implementations, the bounds of an interval are represented by FP numbers, and IA operations use directed roundings to make sure that the exact result remains contained in the output interval.If we use the same precision in conventional FP computation and IA, the result of the FP computation is also included in the interval result: Ŷ ∈ Y.The error Eq. ( 2) is bounded by Unfortunately, this bound cannot be computed because Y is not known.However, we can use the inequalities sup |z| ⩽ 2 rad(z) ⩽ 2 sup |z| , where 2 rad([z, z]) = z − z, which apply when 0 ∈ z.Applying this to each component of Y −Y, which contains 0, yields where In (6), the middle part of the inequality is computable, and it cannot be more than twice the value of the actual error bound.
Therefore the error obtained with IA computations is always larger than the errors obtained through FP computations, but of course knowing in which extent is difficult and highly dependent on the kind of computation being performed.

V. NUMERICAL INVESTIGATIONS
We now test the two approaches presented in sections III and IV by means of intensive numerical computations.The Cooley-Tukey algorithm is implemented in Julia (v1.7.2).Julia's generic programming allows our code to be implemented for all formats consistently, which helps fairer performance comparisons.All tests are performed in binary64/doubleprecision arithmetic ( p = 53) for FP and IA approximations, and we used complex multiplications with FMA, so that ρ = 2u.IA is implemented with the IntervalArithmetic library1 .To estimate errors, we need "exact" values to serve as reference.To that purpose, we used Johansson's Arblib library [15] (through Julia's Arblib library 2 ) to obtain highprecision enclosure of the exact values.

A. Evaluation of errors
Given an input vector X of FP coefficients, we first compute an accurate enclosure of Y using Arblib in high-precision ( p = 256), then we compute a FP (resp.IA) approximation Ŷ (resp.Y ) in double-precision arithmetic ( p = 53).To perform these computations, we obtained beforehand high-precision enclosures of the coefficients ω k N , from which we were able to deduce correctly-rounded FP approximations ωk N , and tight enclosures ω k N .Knowing Y allows the computation of usually inaccessible errors , is also computed.This procedure is applied to N stat = 65 536 randomly-chosen input vectors, and we retain the maximum value of each error e FP n , e IA n , and r IA n .

B. Results
Figure 2 shows the maximum values obtained for e FP n , e IA n , and r IA n over a random set of 65 536 inputs and for sizes 2 n with n = 1, 2, . . ., 13, allowing the following observations: • All computed FP approximations have error smaller than the "bad case", which clearly is "exceptionally bad".for increasing values of the number of samples N stat , up to 2 24 .It still slowly increases with N stat , indicating that larger errors in the FP calculation are still of non-negligible probability.
Fig. 4 shows the timing ratios between FP and IA arithmetic (and between FP and Arb arithmetic) on the environment we used (Julia).The drop for FFTs of size around 2 11 probably comes from cache effects: for these values what we measure is more the delay of memory transfers than the delay of arithmetic operations.This indicates that for very large FFTs, the time penalty due to IA is between a factor of 2 and 10.

VI. CONCLUSION
Our experiments show that, despite the inherent intervalgrowth of interval arithmetic, local error bounds obtained through interval arithmetic are significantly smaller than the state-of-the-art global error bounds.As a consequence, unless either IA is slow on the platform being used or the global bound is smaller than the error threshold required by the considered application, computing FFTs in IA can be a interesting alternative.Moreover, tests with various different interval libraries may allow us to refine our results.
Fig.2.In binary64/double-precision arithmetic: best known bound bn (Eq.4), bad case wn (Eq.5), and maximum errors e FP n , e IA n , and r IA n obtained over a sample of 65 536 randomly-chosen inputs. 1113 2 15 2 17 2 19 2 21 2 23 As expected, the computed IA approximations have much larger errors than the FP approximations, and the obtained local bound is slightly larger.• Nertheless, the local, IA-based, bounds remain smaller than the global bounds.In double-precision, the number of possible inputs makes exhaustive tests impossible, hence a natural question is whether we performed enough tests to get a reasonable estimate of the errors one meets in practice.Figure3shows e FP • Fig. 4. Timing ratios between computation with high-precision Arb (p = 256) / FP arithmetic and IA / FP arithmetic for increasing input sizes.