¿Cómo implementar una versión multidimensional de NymPy convolve?

Después de esta publicación, se me recomendó hacer una pregunta diferente basada en MCVE . Mi objective es implementar la convolución de NumPy para matrices de entrada de forma arbitraria. Tenga en cuenta que estoy usando el término convolución en el contexto del producto de Cauchy de series multivariables de potencia (multiplicación polinomial multivariable). Las funciones de SciPy como signal.convolve , ndimage.convolve o ndimage.filters.convolve no me funcionan como he explicado aquí .

Considere dos matrices numéricas 2D y A no cuadradas A y 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) 

por ejemplo:

 [[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]] 

Ahora puedo calcular los elementos de C=convolve(A,B) con conv(A,B,K) :

 def crop(A,D1,D2): return A[tuple(slice(D1[i], D2[i]) for i in range(A.ndim))] def sumll(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 sumll(np.multiply(crop(A,np.maximum(D0,np.minimum(A.shape,KB.shape)) \ ,np.minimum(A.shape,K)), \ flipall(crop(B,np.maximum(D0,np.minimum(B.shape,KA.shape)) \ ,np.minimum(B.shape,K))))) 

Ejemplo de ejemplo para K=np.array([0,0])+1 , conve(A,B,K) resultados 1*2=2 , para K=np.array([1,0])+1 resultados 5*1+2*6=17 , para K=np.array([0,1])+1 es 2*4+1*2=10 y para K=np.array([1,1])+1 da 4*5+6*2+1*1+1*2=36 :

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

Ahora, si supiera la dimensión de A y B , podría anidar algunos bucles para poblar la C , pero ese no es el caso para mí. ¿Cómo puedo usar la función conv para rellenar el C ndarray con una forma de C.shape=A.shape+B.shape-1 sin usar para bucles?

Existe una multitud de funciones de convolución n-dimensionales en scipy.ndimage y astropy . A ver si podemos usar alguno de ellos.

Primero necesitamos algunos datos para comparar. Así que vamos a expandir el espacio de entrada:

 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]]) 

y calcula tu convolución sobre este espacio:

 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.]]) 

Bien, ahora que sabemos qué valores esperar, veamos si podemos obtener scipy y astropy para darnos los mismos valores:

 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.]]) 

Como puede ver, es solo cuestión de elegir la normalización y el relleno correctos, y simplemente puede usar cualquiera de estas bibliotecas.

Recomiendo usar astropy.convolution.convolve_fft , ya que (al estar basado en FFT) es probablemente el más rápido.