Small value artifacts from applying scipy.signal.convolve

Question

This Python 3.11.5 script:

import numpy as np
from skimage.io import imread
from scipy.signal import convolve

image = np.flipud(imread('conv-test.bmp').astype(np.float32))
con = convolve(image, np.ones((3, 3))/9, mode='valid')
print(image.min())
print(np.logical_and(image > 0, image < 1).any())
print(np.logical_and(con > 0, con < 0.00001).any())

produces:

0.0
False
True

How is it possible to get positive values smaller than 0.00001? Shouldn't 1/9 be the smallest positive value produced by this convolution?

Answer 1

The function scipy.signal.convolve has a method parameter that allows you to select the method that is used to compute the convolution. The available methods are the direct method and the FFT. If you don't explicitly set a method, it will choose a method based on heuristics that in general result in choosing the faster method. In your case, it appears to have selected the FFT method. The FFT is subject to normal floating point imprecision that can result in a value that should theoretically be zero being instead very small but nonzero. Here's a small 1D example:

In [76]: x = np.zeros(18)

In [77]: x[9:] = 1

In [78]: convolve(x, np.ones(3)/3, method='fft', mode='valid')
Out[78]: 
array([-1.49880108e-16, -8.88178420e-17,  0.00000000e+00,  0.00000000e+00,
       -8.88178420e-17, -1.11022302e-16, -8.88178420e-17,  3.33333333e-01,
        6.66666667e-01,  1.00000000e+00,  1.00000000e+00,  1.00000000e+00,
        1.00000000e+00,  1.00000000e+00,  1.00000000e+00,  1.00000000e+00])

Theoretically, the first seven values in that result should all be 0. We can get that result by selecting the direct method:

In [79]: convolve(x, np.ones(3)/3, method='direct', mode='valid')
Out[79]: 
array([0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.33333333, 0.66666667, 1.        ,
       1.        , 1.        , 1.        , 1.        , 1.        ,
       1.        ])