python - 在 numpy 中更有效地矢量化此卷积类型循环

Question

我必须执行以下类型的许多循环

for i in range(len(a)):
    for j in range(i+1):
        c[i] += a[j]*b[i-j]

其中 a 和 b 是短数组(大小相同，大约在 10 到 50 之间)。这可以使用卷积有效地完成:

import numpy as np
np.convolve(a, b) 

但是，这给了我完整的卷积(即，与上面的 for 循环相比，向量太长了)。如果我在卷积中使用“相同”选项，我会得到中心部分，但我想要的是第一部分。当然，我可以从完整向量中去掉不需要的部分，但如果可能的话，我想去掉不必要的计算时间。有人可以建议更好的循环矢量化吗？

Answer 1

您可以在 Cython 中编写一个小的 C 扩展:

# cython: boundscheck=False
cimport numpy as np
import numpy as np  # zeros_like

ctypedef np.float64_t np_t
def convolve_cy_np(np.ndarray[np_t] a not None,
                   np.ndarray[np_t] b not None,
                   np.ndarray[np_t] c=None):
    if c is None:
       c = np.zeros_like(a)
    cdef Py_ssize_t i, j, n = c.shape[0]
    with nogil:
        for i in range(n):
            for j in range(i + 1):
                c[i] += a[j] * b[i - j]
    return c

与我机器上的 np.convolve(a,b)[:len(a)] 相比，它在 n=10..50 上表现良好。

这似乎也是 numba 的工作.