scipy.signal.fftconvolve is not commutative

[lego-engineer]

scipy.signal.fftconvolve is not commutative in 'same' mode with input arrays of different lengths. It mentions this in the documentation (see below), but it seems like swapping the input arrays, warning the user, or erroring would also be appropriate.

scipy.signal.fftconvolve Documentation:

in2 : array_like
        Second input. Should have the same number of dimensions as `in1`;
        if sizes of `in1` and `in2` are not equal then `in1` has to be the
        larger array.

Reproducing code example:

import numpy as np
import scipy.signal as sig

modes = ("valid", "full", "same")
a = np.random.random(50)
b = np.random.random(5)
for mode in modes:
    print(mode)
    assert np.isclose(sig.fftconvolve(a, b, mode), sig.fftconvolve(b, a, mode)).all()
    assert np.isclose(np.convolve(a, b, mode), sig.fftconvolve(a, b, mode)).all()

Error message:

valid
same
Traceback (most recent call last):
  File "<stdin>", line 3, in <module>
  File "/usr/local/lib/python3.5/dist-packages/numpy/core/numeric.py", line 2524, in isclose
    return within_tol(x, y, atol, rtol)
  File "/usr/local/lib/python3.5/dist-packages/numpy/core/numeric.py", line 2510, in within_tol
    return less_equal(abs(x-y), atol + rtol * abs(y))
ValueError: operands could not be broadcast together with shapes (50,) (5,)

Scipy/Numpy/Python version information:

1.2.1 1.16.2 sys.version_info(major=3, minor=5, micro=3, releaselevel='final', serial=0)

[endolith]

but it seems like swapping the input arrays, warning the user, or erroring would also be appropriate.

Is there a good argument for why it should behave in any of these particular ways?

[lego-engineer]

Mathematically, convolution is a commutative operation, so it would make sense for fftconvolve to also be commutative. Alternatively, to maintaining a somewhat uniform interface between numpy.convolve and scipy.signal.convolve.

Anecdotally, I first noticed this working on a SDR project with numpy.convolve and wanted to see if I could get any speed up with scipy.signal.fftconvolve. I briefly skimmed the docs and saw it had roughly the same interface as numpy.convolve so I switched them out and it kept getting garbage output half the time. I quickly learned that fftconvolve was not commutative and added a swap step.

I know my use case was relatively simple, but I don't see how a warning could hurt and, in my case, it would have helped.

[endolith]

Mathematically, convolution is a commutative operation, so it would make sense for fftconvolve to also be commutative.

It is commutative, though:

scipy.signal.fftconvolve([1, 2], [3, 4, 5], 'full')
Out[13]: array([ 3., 10., 13., 10.])

scipy.signal.fftconvolve([3, 4, 5], [1, 2], 'full')
Out[14]: array([ 3., 10., 13., 10.])

But same mode truncates the output so that it matches the shape of the first input:

scipy.signal.fftconvolve([3, 4, 5], [1, 2], 'same')
Out[15]: array([ 3., 10., 13.])

scipy.signal.fftconvolve([1, 2], [3, 4, 5], 'same')
Out[16]: array([10., 13.])

This is typically used for things like image filtering, where you want the output to be the same dimensions as the input, and unshifted (if the kernel is symmetrical). So if the kernel is identity-like, the output should be unchanged/unshifted:

scipy.signal.fftconvolve([1, 2, 3, 4], [0, 1, 0], 'same')
Out[17]: array([1., 2., 3., 4.])

scipy.signal.fftconvolve([1, 2], [0, 1, 0], 'same')
Out[18]: array([1., 2.])

scipy.signal.fftconvolve([1, 2, 3, 4], [0.01, 1, 0.01], 'same')
Out[19]: array([1.02, 2.04, 3.06, 4.03])

It is still commutative if the inputs are both the same shape:

scipy.signal.fftconvolve([1, 2, 3], [4, 5, 6], 'same')
Out[4]: array([13., 28., 27.])

scipy.signal.fftconvolve([4, 5, 6], [1, 2, 3], 'same')
Out[5]: array([13., 28., 27.])

But when the inputs are not the same shape, how would you change this so that the truncated output is commutative but also matches the shape of the first input? That doesn't seem possible to me?

[lego-engineer]

I agree that one could not make it match the first input and have it be commutative.

While I agree that the current definition is convenient some filtering applications, I am more accustomed to "same" indicating a convolution size of max(len(in1), len(in2)). In my experience, this alternate definition works just as well for filtering as signals are usually significantly larger than filter kernels. However, even with the current definition, I think a warning would be useful. People who mean to operate that way would just ignore the warning, people who don't mean to would be warned that they should swap their arguments.

Moreover, I think there is still a good argument for maintaining a more uniform interface between numpy.convolve, scipy.signal.convolve, and scipy.signal.fftconvolve, especially when scipy.signal.convolve is calling scipy.signal.fftconvolve under the hood.

[endolith]

I am more accustomed to "same" indicating a convolution size of max(len(in1), len(in2)).

Just because you're used to it from some other software? Or is there some actual benefit to that method, over just matching the shape of the first input in a predictable way? I don't understand why it would be desirable to change behavior based on the input shape.

Also, which dimension should we be comparing when the input is N-dimensional?

scipy.signal.convolve([[1, 2, 3], [4, 5 ,6]],
                      [[1, 2], [3, 4], [5, 6]], 'same')
                      
Out[6]: 
array([[ 7, 23, 33],
       [17, 47, 65]])

signals are usually significantly larger than filter kernels.

But not always, right?

However, even with the current definition, I think a warning would be useful. People who mean to operate that way would just ignore the warning, people who don't mean to would be warned that they should swap their arguments.

Wouldn't it make more sense to post this on numpy, since their behavior is more likely to be unexpected? Scipy's behavior is consistent with Matlab's conv/conv2/convn/etc. that people are often familiar with.

Moreover, I think there is still a good argument for maintaining a more uniform interface between numpy.convolve, scipy.signal.convolve, and scipy.signal.fftconvolve,

I wouldn't be opposed to adding another mode for compatibility purposes, but I'm not sure how the ND case would be handled.

especially when scipy.signal.convolve is calling scipy.signal.fftconvolve under the hood.

Don't those already behave consistently?

[endolith]

Actually I was wrong, Matlab does always copy the shape of the first input, but the center of the output can be offset by one relative to SciPy:

Matlab example:

u = [-1 2 3 -2 0 1 2];
v = [2 4 -1 1];
w = conv(u,v,'same')

w = 1×7

    15     5    -9     7     6     7    -1

While SciPy outputs array([ 0, 15, 5, -9, 7, 6, 7])

The full output is 

          ------------------------+------- Matlab
          ↓                       ↓
[-2,  0, 15,  5, -9,  7,  6,  7, -1,  2]
      ↑                       ↑
      ------------------------+----------- SciPy

Might be worthwhile to have a mode='same_matlab' or something for compatibility with existing algorithms?

[lego-engineer]

I wouldn't by any means call the numpy method unpredictable but I do see how it deviates from Matlab and I know many people still use Matlab for prototyping. However, I see now that I was mistaken in scipy.fft.convolve's operation. Briefly looking at the source, I saw it was calling numpy's convolve function and thought it had the same output. However, since scipy maintains a convolution consistent convolution interface, I think my previous suggestions are less appropriate.

If there is interest, I could make up a PR for "same_numpy" mode, but it seems like that may be cluttering up the interface without a significant benefit. We also may consider adding a note to the scipy documentation explicitly calling out that "same" numpy and scipy convolution aren't always equivalent. I would argue that since scipy is built on top of numpy, it makes more sense to have such a note in scipy than numpy.

I think you have found an interesting bug in Matlab / Scipy "same" compatibility that likely should be resolved (I think the inclusion of a new mode seems reasonable).

[endolith]

Documentation is definitely a good idea.

I'd probably be supportive of other compatibility modes, too, It's always nice to be able to just use something and trust that someone else already figured out the corner cases and it will just work like it's supposed to. Not sure how to name them though.

[cjblocker]

The difference only appears when the size of the second argument is even, otherwise, MATLAB and SciPy give the same result. In the even case, the center is ambiguous, but it's not in the odd case.

I think this should definitely be in the documentation, as this is the convolution routine most people porting from MATLAB use since it is the only one that offers the 'same' mode. When I read the above comment, I was a little scared because I had ported a lot of code a long time ago from MATLAB, but fortunately, all of my kernels are odd sized.

Assuming the difference is not considered a bug, I don't think a compatibility mode is necessary. If you need the same output, MATLAB compatibility can easily be achieved by post-padding your even filter dimensions with zeros, like so:

import numpy as np
from scipy.signal import convolve as conv

def pad_to_odd(kernel):
    return np.pad(kernel, [(0,1) if dim%2==0 else (0,0) for dim in kernel.shape], mode='constant')

u = np.array([-1, 2, 3, -2, 0, 1, 2, 1])
v = np.array([2, 4, -1, 1])
w = conv(u,pad_to_odd(v),mode='same')
w  
    array([15,  5, -9,  7,  6,  9,  3,  1])

which is what MATLAB gives

u = [-1 2 3 -2 0 1 2 1];
v = [2 4 -1 1 ];
w = conv2(u,v,'same')
w =

    15     5    -9     7     6     9     3     1

[endolith]

@cjblocker

If you need the same output, MATLAB compatibility can easily be achieved by post-padding your even filter dimensions with zeros,

But it's even easier to just use a compatibility mode that handles it for you, and you don't have to worry about making mistakes in your implementation.

[larsoner]

My vote would be just to document the difference rather than adding a new mode

[ilayn]

My vote would be just to document the difference rather than adding a new mode

+1 There is already a lot of confusion even in matlab funcs about this.