How to implement a multidimensional version of NymPy convolve?

Question

Following this post I was advised to ask a different question based on MCVE. My objective is to implement the NumPy's convolve for arbitrary shaped input arrays. Please consider that I'm using the term convolution in the context of Cauchy product of multivariate power series (multivariable polynomial multiplication). SciPy functions such as signal.convolve, ndimage.convolve or ndimage.filters.convolve do not work for me as I have explained here.

Consider two non-square 2D NumPy arrays A and B:

D1=np.array([4,5])
D2=np.array([2,3])
A=np.random.randint(10,size=D1)
B=np.random.randint(10,size=D2)

for example:

[[1 4 4 2 7]
 [6 1 7 5 3]
 [1 4 3 4 8]
 [7 5 8 3 3]]
[[2 2 3]
 [5 2 9]]

Now I'm able to calculate the elements of the C=convolve(A,B) with conv(A,B,K):

def crop(A,D1,D2):
    return A[tuple(slice(D1[i], D2[i]) for i in range(A.ndim))]

def sumall(A):
    sum1=A
    for k in range(A.ndim):
        sum1 = np.sum(sum1,axis=0)
    return sum1

def flipall(A):
    return A[[slice(None, None, -1)] * A.ndim]

def conv(A,B,K):
    D0=np.zeros(K.shape,dtype=K.dtype)
    return sumall(np.multiply(crop(A,np.maximum(D0,np.minimum(A.shape,K-B.shape)) \
                           ,np.minimum(A.shape,K)), \
                      flipall(crop(B,np.maximum(D0,np.minimum(B.shape,K-A.shape)) \
                           ,np.minimum(B.shape,K)))))

Fow example for K=np.array([0,0])+1, conve(A,B,K) results 1*2=2, for K=np.array([1,0])+1 results 5*1+2*6=17, for K=np.array([0,1])+1 is 2*4+1*2=10 and for K=np.array([1,1])+1 gives 4*5+6*2+1*1+1*2=36:

[[2 10 ...]
 [17 36 ...]
...  ]]

now if I knew the dimension of A and B I could nest some for loops to populate the C, but that's not the case for me. How can I use the conv function to populate the C ndarray with a shape of C.shape=A.shape+B.shape-1 without using for loops?

There must be some typo in your expected return values. Can you please post some actual code and some actual returned values? So that I can just copy and paste it into my editor. — Nils Werner, Aug 13 '18 at 11:55
@NilsWerner They are indeed the actual results results I copy pasted from my notebook. please see [this gist](https://gist.github.com/Foadsf/055ef75dd0109e47604b800c90ad5326). — Foad S. Farimani, Aug 13 '18 at 12:00
I am pretty sure `2*2 != 2`, and don't post Gists but amend your post. — Nils Werner, Aug 13 '18 at 12:00
dear @NilsWerner the gist is just a replica of the post above, nothing else has been mentioned there. — Foad S. Farimani, Aug 13 '18 at 12:01

Nils Werner · Answer 1 · 2018-08-20T09:21:00.627

2

There exist a multitude of n-dimensional convolution functions in scipy.ndimage and astropy. Let's see if we can use any of them.

First we need some data to compare against. So let's span up the input space:

d0, d1 = np.array(A.shape) + np.array(B.shape) - 1
input_space = np.array(np.meshgrid(np.arange(d0), np.arange(d1))).T.reshape(-1, 2)
# array([[0, 0],
#        [0, 1],
#        [0, 2],
#        [0, 3],
#        [0, 4],
#        [0, 5],
#        [0, 6],
#        [1, 0],
#        [1, 1],
#        ...
#        [4, 5],
#        [4, 6]])

and calculate your convolution over this space:

out = np.zeros((d0, d1))
for K in input_space:
    out[tuple(K)] = conv(A, B, K + 1)

out
# array([[  2.,  10.,  19.,  24.,  30.,  20.,  21.],
#        [ 17.,  36.,  71.,  81., 112.,  53.,  72.],
#        [ 32.,  27., 108.,  74., 121.,  79.,  51.],
#        [ 19.,  46.,  79.,  99., 111.,  67.,  81.],
#        [ 35.,  39., 113.,  76.,  93.,  33.,  27.]])

Okay, now that we know what values to expect, lets see if we can get scipy and astropy to give us the same values:

import scipy.signal
scipy.signal.convolve2d(A, B)  # only 2D!
# array([[  2.,  10.,  19.,  24.,  30.,  20.,  21.],
#        [ 17.,  36.,  71.,  81., 112.,  53.,  72.],
#        [ 32.,  27., 108.,  74., 121.,  79.,  51.],
#        [ 19.,  46.,  79.,  99., 111.,  67.,  81.],
#        [ 35.,  39., 113.,  76.,  93.,  33.,  27.]])

import astropy.convolution
astropy.convolution.convolve_fft(
    np.pad(A, pad_width=((1, 0), (1, 1)), mode='constant'),
    B,
    normalize_kernel=False
)
# array([[  2.,  10.,  19.,  24.,  30.,  20.,  21.],
#        [ 17.,  36.,  71.,  81., 112.,  53.,  72.],
#        [ 32.,  27., 108.,  74., 121.,  79.,  51.],
#        [ 19.,  46.,  79.,  99., 111.,  67.,  81.],
#        [ 35.,  39., 113.,  76.,  93.,  33.,  27.]])
astropy.convolution.convolve(
    np.pad(A, pad_width=((1, 0), (1, 1)), mode='constant'),
    np.pad(B, pad_width=((0, 1), (0, 0)), mode='constant'),
    normalize_kernel=False
)
# array([[  2.,  10.,  19.,  24.,  30.,  20.,  21.],
#        [ 17.,  36.,  71.,  81., 112.,  53.,  72.],
#        [ 32.,  27., 108.,  74., 121.,  79.,  51.],
#        [ 19.,  46.,  79.,  99., 111.,  67.,  81.],
#        [ 35.,  39., 113.,  76.,  93.,  33.,  27.]])

import scipy
scipy.ndimage.filters.convolve(
    np.pad(A, pad_width=((0, 1), (0, 2)), mode='constant'),
    B,
    mode='constant',
    cval=0.0,
    origin=-1
)
# array([[  2.,  10.,  19.,  24.,  30.,  20.,  21.],
#        [ 17.,  36.,  71.,  81., 112.,  53.,  72.],
#        [ 32.,  27., 108.,  74., 121.,  79.,  51.],
#        [ 19.,  46.,  79.,  99., 111.,  67.,  81.],
#        [ 35.,  39., 113.,  76.,  93.,  33.,  27.]])
scipy.ndimage.filters.convolve(
    np.pad(A, pad_width=((1, 0), (1, 1)), mode='constant'),
    B,
    mode='constant',
    cval=0.0
)
# array([[  2.,  10.,  19.,  24.,  30.,  20.,  21.],
#        [ 17.,  36.,  71.,  81., 112.,  53.,  72.],
#        [ 32.,  27., 108.,  74., 121.,  79.,  51.],
#        [ 19.,  46.,  79.,  99., 111.,  67.,  81.],
#        [ 35.,  39., 113.,  76.,  93.,  33.,  27.]])

As you can see, it is just a matter of chosing the right normalization and padding, and you can simply use any of these libraries.

I recommend using astropy.convolution.convolve_fft, as it (being FFT-based) is probably the fastest.

edited Aug 20 '18 at 09:21

answered Aug 13 '18 at 12:56

Nils Werner

34,832
7
76
98

this is awesome. You answered my original question. I'm gonna try them and come here. I appreciate your support. – Foad S. Farimani Aug 13 '18 at 13:03
I think now I understand what you have done. functions above are actually meant for image processing so you you have expanded them with `numpy.pad` and then used the scipy or astropy functions. but this requires the shape of `A` and `B` to be known. if we automate that part then I get what I need. – Foad S. Farimani Aug 13 '18 at 13:25
How can the shapes of `A` and `B` ever not be known? – Nils Werner Aug 13 '18 at 13:28
My final goal is to multiply multivariate finite power series. at some point I have to increase the degree of the serie for specific variable. Although doing so manually is possible it will be inefficient. It would be ideal to construct the Convolve function to take `A` and `B` and gives the result. – Foad S. Farimani Aug 13 '18 at 13:33
1

Thats what all of these libraries do! You give them `A` and `B` (maybe pad some zeros here and there), and you get the result. – Nils Werner Aug 13 '18 at 13:33
1

Did this answer help you to solve your problem? If so, please remember to [mark it as accepted](https://meta.stackexchange.com/questions/5234/how-does-accepting-an-answer-work). – Nils Werner Aug 20 '18 at 09:00
You answer indeed helped me to learn that I was wrong about the existing functions and how one may modify them to calculate the convolution. For that matter I upvoted it. But because it doesn't answer the question of how to populate the `C` function using the `conv` function this can't be an accepted answer. unless I change the question in a way to fit this answer? – Foad S. Farimani Aug 20 '18 at 09:06
Just do `C = astropy.convolution.convolve_fft(...)`. If that's not what you want you should explicitly state some example code that assigns values to `C` the way you want! – Nils Werner Aug 20 '18 at 09:09
:) yes of course. The title of the question is `How to implement a multidimensional version of NymPy convolve?`, What I have had in mind is to have the same function as [numpy convolve](https://docs.scipy.org/doc/numpy/reference/generated/numpy.convolve.html), with panning and everything else automatically done, as mentioned above. – Foad S. Farimani Aug 20 '18 at 09:12
All of the aforementioned convolve functions work on n dimensions, with the exception of `scipy.signal.convolve2d`! What more do you want? – Nils Werner Aug 20 '18 at 09:13
I want a convolve function, just like [NumPy convolve](https://docs.scipy.org/doc/numpy/reference/generated/numpy.convolve.html) which takes two ndarrays and spits out the result, sparing the user from panning and other tidius considerations. – Foad S. Farimani Aug 20 '18 at 09:17
Have you even tried any of the mentioned functions? None of them require the user to do any panning! – Nils Werner Aug 20 '18 at 09:18
sorry bad writing. I meant padding. Yes indeed I tried them the day I wrote the post. – Foad S. Farimani Aug 20 '18 at 09:21
You can write 100 lines to create a slow and untested n-dimensional convolution. Or you can write 3 lines to adjust your data a little bit so it works with `scipy` or `astropy`, both of which are fast and are known to work. Why do you think the former is the better solution? – Nils Werner Aug 20 '18 at 09:23
It is not my goal to rewrite those functions. Indeed if we could take care of the padding automatically then the problem was solved. It is indeed possible and I'm thinking of implementing it and adding it to your solution. However I was busy with some other tasks... – Foad S. Farimani Aug 20 '18 at 09:25
FYI [here](https://mail.python.org/pipermail/astropy/2018-August/004445.html) I was told that `scipy.ndimage.correlate` with the right `mode` and `cval` parameters will do the job. I haven't checked it yet though. – Foad S. Farimani Aug 20 '18 at 09:31

How to implement a multidimensional version of NymPy convolve?

1 Answers1

Linked