# Fast & accurate convolution algorithm (like FFT) for high dynamic range?

It seems that FFT-based convolution suffers from limited floating-point resolution due to evaluating everything around the roots of unity, as you can see in the $10^{14}$-factor error in this Python code:

from scipy.signal import convolve, fftconvolve
a = [1.0, 1E-15]
b = [1.0, 1E-15]
convolve(a, b)     # [  1.00000000e+00,   2.00000000e-15,   1.00000000e-30]
fftconvolve(a, b)  # [  1.00000000e+00,   2.11022302e-15,   1.10223025e-16]


Are there any fast convolution algorithms that do not suffer from this problem?
Or is direct (quadratic-time) convolution the only way to get an accurate solution?

(Whether such small numbers are significant enough not to chop off is beside my point.)

• Note that convolve() just calls fftconvolve() now, if the input size is large. Specify method='direct' if you want direct. Jun 16 '17 at 21:08
• @endolith: Good point! I just learned that recently but forgot about it here. Jun 17 '17 at 0:01

Disclaimer: I know this topic is older, but if one is looking for "fast accurate convolution high dynamic range" or similar this is one of the first of only a few decent results. I wanna share my insights I got on this topic so it might help somebody in the future. I apologize if I might use the wrong terms in my answer, but everything I found on this topic is rather vague and lead to confusing even in this thread. I hope the reader will understand anyway.

Direct convolution is mostly accurate to machine precision for each point, i.e. the relative error is usually roughly or close to 1.e-16 for double precision for each point of the result. Each point has 16 correct digits. Rounding errors can significant for untypically large convolutions though, and strictly speaking one should be careful with cancellation and use something like Kahan summation and high enough precision data types, but in practice the error as almost always optimal.

The error of a FFT convolution apart from rounding errors is "global relative" error, meaning the error in each point depends on the machine precision and the peak value of the result. For example if the peak value of the result is 2.e9, then the absolute error in each point is $$2\cdot10^9\cdot10^{-16} = 2\cdot10^{-7}$$. So if a value in the result is supposed to be very small, let's say $$10^{-9}$$, the relative error in that point can be huge. FFT convolution is basically useless if you need small relative errors in the tail of your result, e.g. you have a somewhat exponential decay of your data and need accurate values in the tail. Interestingly if FFT convolution is not limited by that error, it has much smaller rounding errors compared to direct convolution, since you obviously do less additions/multiplications. This is actually why people often claim FFT convolution is more accurate, and they are almost right in some sense, so they can be quite adamant.

Unforunately there is no easy universal fix to get fast and accurate convolutions, but depending on your problem there might be one... I have found two:

If you have smooth kernels which can be approximated well by a polynomial in the tail, then the black-box Fast Multipole Method with Chebyshev interpolation might be interesting for you. If your kernel is "nice" this works actually perfectly: you get both linear (!) computational complexity and machine precision accuracy. If this fits your problem you should use it. It's not easy to implement however.

For some specific kernels (convex functions I think, usually from probability densities) you can use an "exponential shift" to get optimal error in some part of the tail of the result. There is a PHD thesis and a github with a python implementation using that systematically, and the author cals it accurate FFT convolution. In most cases this is not super useful however, since either it regresses back to direct convolution or you can use FFT convolution anyway. Although the code does it automatically, which is nice of course.

--------------------EDIT:--------------------

I looked a little bit at the Karatsuba algorithm (I actually made a small implementation), and to me it looks like it has usually similar error behavior like the FFT convolution, i.e. you get an error relative to the peak value of the result. Due to the divide and conquer nature of the algorithm some values in the tail of the result actually have better error, but I don't see an easy systematic way to tell which ones or in any case how to use this observation. Too bad, at first I thought Karatsuba might be something useful in-between direct and FFT convolution. But I don't see common use cases where Karatsuba should be prefered over the common two convolution algorithms.

And to add to the exponential shift I mentioned above: There are many cases where you can use it to improve the result of a convolution, but again it's not a universal fix. I actually use this together with FFT convolution to get pretty good results (in the general case for all inputs: at the worst same error as normal FFT convolution, at best relative error in each point to machine precision). But again, this only really works nicely for specific kernels and data, but for me both kernel and data or somewhat exponential in decay.

• +1 welcome & thank you very much for posting this! :) Oct 24 '18 at 9:40
• Wow! i learned something also and that is a new term for something i've been doing since 1993. this Kahan summation algorithm appears to be exactly the same as what i had been calling noise shaping with a zero in the noise-to-output transfer function placed right at DC or the zero is placed at $z = 1$ on the $z$ plane. Randy Yates called it "fraction saving", which is a concise generic name for it. i wonder who mr/ms Kahan is and when this is credited. Oct 25 '18 at 20:34
• The original publication of Kahan seems to be from 1964.
– oli
Oct 25 '18 at 20:50
• it's today's surprise. Actually a while a go @DanBoschen had asked a dsp puzzle, considering the dynamic range of floating point numbers, which actually was about this same concept of adding very small numbers to very large numbers... Oct 25 '18 at 21:25

One candidate is the Karatsuba algorithm, which runs in $O\big(N^{\log_23}\big) \approx O\big(N^{1.5849625}\big)$ time. It's not transform-based. There's also some code with an explanation in the Music-DSP Source Code Archive, which looks like an independent discovery of a similar algorithm.

Testing a Python implementation of the Karatsuba algorithm (installed by sudo pip install karatsuba) using the numbers in your question shows that even with 64-bit floating point numbers the relative error is large for one of the output values:

import numpy as np
from karatsuba import *
k = make_plan(range(2), range(2))
l = [np.float64(1), np.float64(1E-15)]
np.set_printoptions(formatter={'float': lambda x: format(x, '.17E')})
print "Karatsuba:"
print(k(l, l)[0:3])
print "Direct:"
print(np.convolve(l, l)[0:3])


which prints:

Karatsuba:
[1.0, 1.9984014443252818e-15, 1.0000000000000001e-30]
Direct:
[1.00000000000000000E+00 2.00000000000000016E-15 1.00000000000000008E-30]

• There is a extra ] at the end of the link to Karatsuba algorithm
– user1932
Jun 17 '17 at 9:03
• +1 because it's brilliant and it had never occurred to me that Karatsuba was a convolution algorithm, but it would be nice if you could explain why it should solve this problem. I can easily see it for the 2x2 case, but in the general recursive setting I don't see why it should fix this issue. It would seem plausible to me that it might not even be fixable in general, but I don't know. Jun 17 '17 at 9:38
• @OlliNiemitalo: Well the easy way to explain it is that I want the relative error to be low compared to direct $O(n^2)$ convolution. (Any reasonable definition of "low" would work here... the relative error I'm getting with FFT is like $10^{14}$ which isn't low by any definition.) Jun 17 '17 at 10:10
• IEEE doubles only have a precision of 15 to 16 decimal digits in the general case. So 1e-14 is a reasonable size error for a sequence of some number of arithmetic operations (unless you pick a few magic values). Jun 17 '17 at 21:59
• If you've ever designed a floating point adder, you will know that the exponent is determined by the mantissa result during normalization. You picked numbers that produce an unlikely narrow mantissa. Jun 18 '17 at 3:37

Rather than scrapping the fast convolution algorithm, why not use an FFT with a higher dynamic range?

An answer to this question shows how to use the Eigen FFT library with boost multiprecision.

I believe that the Cordic algorithm's precision can be extended as far as you want, if so use an integer DFT and a word length appropriate to your problem.

The same would be true with direct convolution, use very long integers.

Quadratic time convolution to get a DFT result is usually less accurate (can incur more finite quantization numerical noise, due to a deeper layering of arithmetic steps) than the typical FFT algorithm when using the same arithmetic types and operation units.

You might want to try higher precision data types (quad precision or bignum arithmetic).

• Er, this is using the same arithmetic types and operation units isn't it? Clearly it's more accurate. I think the type of noise you're talking about is not the same as the kind I'm talking about. The roots of unity have a magnitude of 1 which means they simply can't represent very small values. This seems not totally related to the question of how noise propagates through the system. Jun 16 '17 at 23:48
• It only seems more accurate in your example because you picked a length and values where the rounding happened to work in your favor. Try a range of much longer convolutions with many more non-zero coefficients with a distribution containing a wide order of magnitudes. Jun 17 '17 at 19:54
• The problem I'm trying to solve has nothing to do with rounding though. That's a different issue I'm not trying to solve. The original examples I had were exactly like what you just said, and they worked just fine with direct convolution but got destroyed by FFT. Jun 17 '17 at 21:18
• Rounding (or other quantization methods) is involved in all finite-precision arithmetic. Some computational results change when rounded, other don't or change less. Jun 17 '17 at 21:34
• I never claimed otherwise. What I just told you was the problem I'm trying to solve has nothing to do with rounding. It's a different problem. I don't care to avoid rounding, but I do care to avoid this problem. Jun 17 '17 at 21:49