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

Question

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.)

Answer 1

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.

Answer 2

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]

Answer 3

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.

Answer 4

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.

Answer 5

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).