python - numpy中矩形子数组的快速顺序添加

Question

我遇到了一个问题，即以向量化形式重写一段代码。下面显示的代码是初始问题的简化说明

K = 20
h, w = 15, 20 
H, W = 1000-h, 2000-w 
q = np.random.randint(0, 20, size=(H, W, K))  # random just for illustration 
Q = np.zeros((H+h, W+w, K))
for n in range(H):
    for m in range(W):
        Q[n:n+h, m:m+w, :] += q[n, m, :]

这段代码需要很长时间才能执行，在我看来允许矢量化实现相当简单。

我知道 numpy 的 s_允许构造切片的函数，而切片又可以帮助代码矢量化。但是因为 Q 中的每个元素都是 q 元素的多次后续添加的结果，所以我发现很难以这种简单的方式进行。

我猜 np.add.at可能有助于处理顺序添加。但我花了很多时间试图让这两个功能为我工作，并决定寻求帮助，因为我经常得到一个

IndexError: failed to coerce slice entry of type numpy.ndarray to integer

对于我所做的任何尝试。

也许还有一些我不知道的 numpy 的魔法，它可以帮助我完成我的任务，但似乎很难用谷歌搜索它。

Answer 1

好吧，您基本上是沿着第一轴和第二轴对滑动窗口求和，这在信号处理领域被称为 convolution .对于 2D 卷积 的两个轴。现在，Scipy 将其实现为 convolve2d并且可以用于沿第三轴的每个切片。

因此，我们将有一个实现，就像这样 -

from scipy.signal import convolve2d

kernel = np.ones((h,w),dtype=int)
m,n,r = q.shape[0]+h-1, q.shape[1]+w-1, q.shape[2]
out = np.empty((m,n,r),dtype=q.dtype)
for i in range(r):
    out[...,i] = convolve2d(q[...,i],kernel)

事实证明，我们可以使用 fftconvolve来自允许我们使用高维数组的同一个 repo。这将使我们以完全矢量化的方式获得输出，就像这样 -

from scipy.signal import fftconvolve

out = fftconvolve(q,np.ones((h,w,1),dtype=int))

运行时测试

函数定义-

def original_app(q,K,h,w,H,W):
    Q = np.zeros((H+h-1, W+w-1, K))
    for n in range(H):
        for m in range(W):
            Q[n:n+h, m:m+w, :] += q[n, m, :]
    return Q        

def convolve2d_app(q,K,h,w,H,W):
    kernel = np.ones((h,w),dtype=int)
    m,n,r = q.shape[0]+h-1, q.shape[1]+w-1, q.shape[2]
    out = np.empty((m,n,r),dtype=q.dtype)
    for i in range(r):
        out[...,i] = convolve2d(q[...,i],kernel)
    return out

def fftconvolve_app(q,K,h,w,H,W):
    return fftconvolve(q,np.ones((h,w,1),dtype=int))

时间和验证-

In [128]: # Setup inputs
     ...: K = 20
     ...: h, w = 15, 20 
     ...: H, W = 200-h, 400-w
     ...: q = np.random.randint(0, 20, size=(H, W, K))
     ...: 

In [129]: %timeit original_app(q,K,h,w,H,W)
1 loops, best of 3: 2.05 s per loop

In [130]: %timeit convolve2d_app(q,K,h,w,H,W)
1 loops, best of 3: 2.05 s per loop

In [131]: %timeit fftconvolve_app(q,K,h,w,H,W)
1 loops, best of 3: 233 ms per loop

In [132]: np.allclose(original_app(q,K,h,w,H,W),convolve2d_app(q,K,h,w,H,W))
Out[132]: True

In [133]: np.allclose(original_app(q,K,h,w,H,W),fftconvolve_app(q,K,h,w,H,W))
Out[133]: True

因此，看起来基于 fftconvolve 的方法在那里做得非常好!