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.
convolve()
just callsfftconvolve()
now, if the input size is large. Specifymethod='direct'
if you want direct. $\endgroup$ – endolith Jun 16 '17 at 21:08