沿轴进行2个二维数组的Numpy卷积

Question

我有两个二维数组。我想沿着轴1进行卷积。np.convolve没有提供axis参数。答案在这里，使用np.apply_along_axis将一个二维数组和一个一维数组进行卷积。但是它不能直接应用于我的情况。问题在这里没有答案。

MWE如下。

import numpy as np

a = np.random.randint(0, 5, (2, 5))
"""
a=
array([[4, 2, 0, 4, 3],
       [2, 2, 2, 3, 1]])
"""
b = np.random.randint(0, 5, (2, 2))
"""
b=
array([[4, 3],
       [4, 0]])
"""

# What I want
c = np.convolve(a, b, axis=1)  # axis is not supported as an argument
"""
c=
array([[16, 20,  6, 16, 24,  9],
       [ 8,  8,  8, 12,  4,  0]])
"""

我知道可以使用np.fft.fft来实现，但似乎这是一个不必要的步骤来完成一件简单的事情。有没有简单的方法来做到这一点？谢谢。

- learner

Answer 1

为什么不直接使用带有zip的列表推导式？

>>> np.array([np.convolve(x, y) for x, y in zip(a, b)])
array([[16, 20,  6, 16, 24,  9],
       [ 8,  8,  8, 12,  4,  0]])
>>> 

或者使用scipy.signal.convolve2d函数：

>>> from scipy.signal import convolve2d
>>> convolve2d(a, b)[[0, 2]]
array([[16, 20,  6, 16, 24,  9],
       [ 8,  8,  8, 12,  4,  0]])
>>> 

- U13-Forward

---

使用zip的列表推导式在存在三维数组且需要一维卷积时无法工作。需要使用两个循环。 convolve2d也存在类似的问题。您正在为OP给出的示例使用一些技巧，但我认为这是一个有用的问题，并且通用答案对社区会更有益处。 - Nagabhushan S N

Answer 2

循环遍历正交维度是否会被认为太丑陋？除非主维度非常短，否则这不会增加太多开销。提前创建输出数组可以确保不需要复制任何内存。

def convolvesecond(a, b):
    N1, L1 = a.shape
    N2, L2 = b.shape
    if N1 != N2:
       raise ValueError("Not compatible")
    c = np.zeros((N1, L1 + L2 - 1), dtype=a.dtype)
    for n in range(N1):
        c[n,:] = np.convolve(a[n,:], b[n,:], 'full')
    return c

对于通用情况（在一对多维数组的第k个轴上进行卷积），我会使用一对辅助函数，将多维问题转换为基本的二维情况：

def semiflatten(x, d=0):
    '''SEMIFLATTEN - Permute and reshape an array to convenient matrix form
    y, s = SEMIFLATTEN(x, d) permutes and reshapes the arbitrary array X so 
    that input dimension D (default: 0) becomes the second dimension of the 
    output, and all other dimensions (if any) are combined into the first 
    dimension of the output. The output is always 2-D, even if the input is
    only 1-D.
    If D<0, dimensions are counted from the end.
    Return value S can be used to invert the operation using SEMIUNFLATTEN.
    This is useful to facilitate looping over arrays with unknown shape.'''
    x = np.array(x)
    shp = x.shape
    ndims = x.ndim
    if d<0:
        d = ndims + d
    perm = list(range(ndims))
    perm.pop(d)
    perm.append(d)
    y = np.transpose(x, perm)
    # Y has the original D-th axis last, preceded by the other axes, in order
    rest = np.array(shp, int)[perm[:-1]]
    y = np.reshape(y, [np.prod(rest), y.shape[-1]])
    return y, (d, rest)

def semiunflatten(y, s):
    '''SEMIUNFLATTEN - Reverse the operation of SEMIFLATTEN
    x = SEMIUNFLATTEN(y, s), where Y, S are as returned from SEMIFLATTEN,
    reverses the reshaping and permutation.'''
    d, rest = s
    x = np.reshape(y, np.append(rest, y.shape[-1]))
    perm = list(range(x.ndim))
    perm.pop()
    perm.insert(d, x.ndim-1)
    x = np.transpose(x, perm)
    return x

（请注意，reshape 和 transpose 不会创建副本，因此这些函数非常快速。）

有了这些函数，通用形式可以写成：

def convolvealong(a, b, axis=-1):
   a, S1 = semiflatten(a, axis)
   b, S2 = semiflatten(b, axis)
   c = convolvesecond(a, b)
   return semiunflatten(c, S1)

- Daniel Wagenaar

Answer 3

一种可能的方法是手动前往傅里叶频谱，然后返回：

n = np.max([a.shape, b.shape]) + 1
np.abs(np.fft.ifft(np.fft.fft(a, n=n) * np.fft.fft(b, n=n))).astype(int)
# array([[16, 20,  6, 16, 24,  9],
#        [ 8,  8,  8, 12,  4,  0]])

- Nils Werner

---

嗨，我已经提到fft方法有点绕。我想要在时间域中全部完成。 - learner

请记住，根据信号的长度，采用FFT方法可能是最快的。 - Nils Werner

"朴素卷积和快速傅里叶变换之间的复杂度差异是众所周知的。您可以始终对信号进行填充，使其长度成为2的幂。" - Nils Werner

@learner 这个问题有点拐弯抹角，能否具体说明一下？ - jtlz2

@jtlz2，我的意思是将其转换为另一个域，进行乘法运算，然后再带回来似乎是不必要的。我想知道是否可以在不进行域移位的情况下使用axis参数完成所有操作。 - learner