How can I accurately perform FFT-based convolution on np.float128 arrays?

Scipy does not seems to fully support the np.float128 data type since I get the following pretty-clear error when using it on my Linux machine supporting this type properly.

DeprecationWarning: dtype=float128 is not supported by correlate and will raise an error in SciPy 1.17.0. Supported dtypes are: boolean, integer, np.float16,np.float32, np.float64, np.complex64, np.complex128

Still results are different and a bit more accurate. My version of Scipy is the 1.15.1 which is the latest release on pypi at the time of writing.

Technically, pyFFTW is supposed to support it since the documentation explicitly states:

Internally, three precisions of FFT are supported. These correspond to single precision floating point, double precision floating point and long double precision floating point, which correspond to numpy‘s float32, float64 and longdouble dtypes respectively (and the corresponding complex types). The precision is decided by the relevant scheme, which is specified by the dtype of the input array.

---

Actual issue

The thing is this is not a precision issue but a numerical stability issue! Indeed, let's check the output with much smaller double-precision numbers:

PYTHONCopied
a = np.array([1e40, 1, 1e40, 1])
print(convolve(a, a, mode='full', method='direct'))   # Correct
print(convolve(a, a, mode='full', method='fft'))      # Wrong
print()

a = np.array([1e5, 1, 1e5, 1])
print(convolve(a, a, mode='full', method='direct'))   # Correct
print(convolve(a, a, mode='full', method='fft'))      # Correct
print()

a = np.array([1e15, 1, 1e15, 1])
print(convolve(a, a, mode='full', method='direct'))   # Correct
print(convolve(a, a, mode='full', method='fft'))      # Wrong, but... Oh!

Here is the detailed output:

[1.e+80 2.e+40 2.e+80 4.e+40 1.e+80 2.e+40 1.e+00]
[1.e+80 0.e+00 2.e+80 0.e+00 1.e+80 0.e+00 0.e+00]

[1.e+10 2.e+05 2.e+10 4.e+05 1.e+10 2.e+05 1.e+00]
[1.e+10 2.e+05 2.e+10 4.e+05 1.e+10 2.e+05 1.e+00]

[1.e+30 2.e+15 2.e+30 4.e+15 1.e+30 2.e+15 1.e+00]
[1.00000000e+30 1.99520395e+15 2.00000000e+30 3.93577438e+15 1.00000000e+30 1.94544573e+15 0.00000000e+00]

This shows that when the exponent is big, the FFT methods is not correct anymore. In fact, the numerical error is already pretty big with input values like 1e15 despite the double precision data type supports values as big as 1e+308. In fact, the last item is wrong even with input values of 1e8 which is ridiculously small compared to what double precision supports! When the exponent is too big (> ~7.97), some values are clamped to 0 and results are clearly not reliable anymore.

The FFT-based algorithm is numerically inaccurate. This is likely due to a catastrophic cancellation due to huge values being subtracted by small values several times during the algorithm. Increasing the precision helps a bit but the few additional digits available in the 80-bit x86-specific long double data type is far from being enough here (like the 64-bit ones are far from being enough for input values being 1e15).

In fact, if we use the np.float128 data type, we can actually see only a slight improvement:

PYTHONCopied
With 1e8:   [1.e+16 2.e+08 2.e+16 4.e+08 1.e+16 2.e+08 1.e+00]   # Correct
With 1e9:   [1.e+18 2.e+09 2.e+18 4.e+09 1.e+18 2.e+09 1.e+00]   # Correct
With 1e10:  [1.e+20 2.e+10 2.e+20 4.e+10 1.e+20 2.e+10 0.e+00]   # Wrong (last item)

I think the issue arise roughly when:

out_min_abs_value / out_max_abs_value < data_type_epsilon

where out_min_abs_value and out_max_abs_value are respectively the minimum/maximum absolute values of the expected output.

---

Solutions

I do not know if this numerical stability issue can be solved or if this is a fundamental issue of using a FFT for the convolution. That being said, note you can increase a bit further the precision by calling the native (C/Fortran) FFTW implementation supporting quad-precision numbers (on Linux with GCC and possibly Clang too). This should help the push the limit a bit further: number with twice more digits. This additional precision comes at a significant price: a significantly slower computation since quad-precision numbers are not native ones. In fact, 80-bit floating-point (FP) numbers are already pretty slow on x86-64 CPUs (and not even supported on other architecture) because the computations are performed on a x87 FP unit which is rather obsolete since 1~2 decades.

A naive solution would be to re-implement your own FFT with arbitrary-precision arithmetic. For example, you can write a native code using the fast GMP library. If you need a really high precision (much more than in the example), you could even write this in gmpy2 in pure-Python since a big portion of the time should be spent in multiplying huge numbers. For numbers as big as in the example, number of ~2800 bits should be enough (that is ~100 bytes). In this case, >80% of the time of gmpy2 operations are actually overheads, so I advise you to write a native code.

A key to avoid this numerical issue could be to avoid huge variation in input numbers. This is AFAIK called preconditioning (and pretty common). A common practice is to compute the values in the log space, but then this is not a convolution anymore in this case due to the FFT... Thus, I do not know whether a preconditioning is possible in this case.