cython卷积码出乎意料的慢

Question

我需要快速计算一个矩阵，其条目是通过将过滤器与每行的向量进行卷积得到的，对结果向量的条目进行二次采样，然后将结果与另一个向量进行点积。具体来说，我想计算

M = [conv(e_j, f)*P_i*v_i ]_{i,j},

其中 i 从 1 到 n 变化，j 从 1 到 m 变化。这里 e_j 是大小为 n 的指示符（行）向量，仅在 j 列中具有一个，f 是长度为 s 的过滤器，P_i 是一个 (n+s-1)-by-k 矩阵，它从卷积，v_i 是长度为 k 的列向量。

计算 M 的每个条目需要 O(n*s) 次操作，因此计算 M 总共需要 O(n*s*n*m)。对于 n=6，m=7，s=3，一个核心我的计算机（8GLOPs）应该能够在大约 0.094 微秒内计算 M。然而，我遵循example given in the Cython documentation 的非常简单的 cython 实现需要超过 2 毫秒来计算具有这些参数的示例。这大约是 4 个数量级的差异！

这是一个包含 Cython 实现和测试代码的 shar 文件。将其复制并粘贴到文件中并在干净的目录中运行“bash ”以获取代码，然后运行“bash ./test.sh”以查看糟糕的性能。

cat > fastcalcM.pyx <<'EOF'

import numpy as np
cimport numpy as np
cimport cython
from scipy.signal import convolve

DTYPE=np.float32
ctypedef np.float32_t DTYPE_t

@cython.boundscheck(False)
def calcM(np.ndarray[DTYPE_t, ndim=1, negative_indices=False] filtertaps, int
        n, int m, np.ndarray[np.int_t, ndim=2, negative_indices=False]
        keep_indices, np.ndarray[DTYPE_t, ndim=2, negative_indices=False] V):
    """ Computes a numrows-by-k matrix M whose entries satisfy
        M_{i,k} = [conv(e_j, f)^T * P_i * v_i],
        where v_i^T is the i-th row of V, and P_i samples the entries from
        conv(e_j, f)^T indicated by the ith row of the keep_indices matrix """

    cdef int k = keep_indices.shape[1]

    cdef np.ndarray M = np.zeros((n, m), dtype=DTYPE)
    cdef np.ndarray ej = np.zeros((m,), dtype=DTYPE)
    cdef np.ndarray convolution
    cdef int rowidx, colidx, kidx

    for rowidx in range(n):
        for colidx in range(m):
            ej[colidx] = 1
            convolution = convolve(ej, filtertaps, mode='full')
            for kidx in range(k):
                M[rowidx, colidx] += convolution[keep_indices[rowidx, kidx]] * V[rowidx, kidx]
            ej[colidx] = 0

    return M

EOF
#-----------------------------------------------------------------------------
cat > test_calcM.py << 'EOF'

import numpy as np
from fastcalcM import calcM

filtertaps = np.array([-1, 2, -1]).astype(np.float32)
n, m = 6, 7
keep_indices = np.array([[1, 3], 
                         [4, 5],
                         [2, 2], 
                         [5, 5], 
                         [3, 4], 
                         [4, 5]]).astype(np.int)
V = np.random.random_integers(-5, 5, size=(6, 2)).astype(np.float32)

print calcM(filtertaps, n, m, keep_indices, V)

EOF
#-----------------------------------------------------------------------------
cat > test.sh << 'EOF'

python setup.py build_ext --inplace
echo -e "%run test_calcM\n%timeit calcM(filtertaps, n, m, keep_indices, V)" > script.ipy
ipython script.ipy

EOF
#-----------------------------------------------------------------------------
cat > setup.py << 'EOF'

from distutils.core import setup
from Cython.Build import cythonize
import numpy

setup(
    name="Fast convolutions",
    include_dirs = [numpy.get_include()],
    ext_modules = cythonize("fastcalcM.pyx")
)

EOF

我认为对 scipy 的 convolve 的调用可能是罪魁祸首（我不确定 cython 和 scipy 可以很好地配合使用），所以我实现了自己的卷积代码，在 Cython 文档中使用了相同的示例，但这导致了整体代码大约慢了 10 倍。

关于如何接近理论上可能的速度的任何想法，或差异如此之大的原因？

Answer 1

一方面，M、eg 和 convolution 的输入不允许快速索引。实际上，您所做的输入并没有什么特别的帮助。

但这没关系，因为你有两个开销。首先是在 Cython 和 Python 类型之间进行转换。如果你想经常将它们传递给 Python，你应该保留无类型数组，以防止需要转换。出于这个原因，只需将其移至 Python 即可加快速度（1ms → 0.65μs）。

然后我分析了它：

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    15                                           def calcM(filtertaps, n, m, keep_indices, V):
    16      4111         3615      0.9      0.1      k = keep_indices.shape[1]
    17      4111         8024      2.0      0.1      M = np.zeros((n, m), dtype=np.float32)
    18      4111         6090      1.5      0.1      ej = np.zeros((m,), dtype=np.float32)
    19                                           
    20     28777        18690      0.6      0.3      for rowidx in range(n):
    21    197328       123284      0.6      2.2          for colidx in range(m):
    22    172662       112348      0.7      2.0              ej[colidx] = 1
    23    172662      4076225     23.6     73.6              convolution = convolve(ej, filtertaps, mode='full')
    24    517986       395513      0.8      7.1              for kidx in range(k):
    25    345324       668309      1.9     12.1                  M[rowidx, colidx] += convolution[keep_indices[rowidx, kidx]] * V[rowidx, kidx]
    26    172662       120271      0.7      2.2              ej[colidx] = 0
    27                                           
    28      4111         2374      0.6      0.0      return M

在您考虑任何事情之前，请先处理convolve。

为什么convolve 很慢？嗯，它有很多开销。 numpy/scipy 通常可以；它最适合大型数据集。如果您知道数组的大小将保持较小，只需在 Cython 中重新实现 convolve。

哦，尝试使用缓冲区语法。将DTYPE[:, :] 用于二维数组，DTYPE[:] 用于一维数组等。这是 memoryview 协议，而且效果更好。在某些情况下，它的开销更大，但这些通常可以解决，而且在大多数其他方面都更好。

---

编辑：

您可以尝试（递归地）内联scipy 函数：

import numpy as np
from scipy.signal.sigtools import _correlateND

def calcM(filtertaps, n, m, keep_indices, V):
    k = keep_indices.shape[1]
    M = np.zeros((n, m), dtype=np.float32)
    ej = np.zeros((m,), dtype=np.float32)

    slice_obj = [slice(None, None, -1)] * len(filtertaps.shape)
    sliced_filtertaps_view = filtertaps[slice_obj]

    ps = ej.shape[0] + sliced_filtertaps_view.shape[0] - 1
    in1zpadded = np.zeros(ps, ej.dtype)
    out = np.empty(ps, ej.dtype)

    for rowidx in range(n):
        for colidx in range(m):
            in1zpadded[colidx] = 1

            convolution = _correlateND(in1zpadded, sliced_filtertaps_view, out, 2)

            for kidx in range(k):
                M[rowidx, colidx] += convolution[keep_indices[rowidx, kidx]] * V[rowidx, kidx]

            in1zpadded[colidx] = 0

    return M

请注意，这使用私有实现细节。

这是为特定维度量身定制的，所以我不知道它是否适用于您的实际数据。但它消除了绝大多数开销。然后，您可以通过再次输入来改进：

import numpy as np
cimport numpy as np
from scipy.signal.sigtools import _correlateND

DTYPE=np.float32
ctypedef np.float32_t DTYPE_t

def calcM(filtertaps, int n, int m, np.int_t[:, :] t_keep_indices, DTYPE_t[:, :] t_V):
    cdef int rowidx, colidx, kidx, k
    cdef DTYPE_t[:, :] t_M
    cdef DTYPE_t[:] t_in1zpadded, t_convolution

    k = t_keep_indices.shape[1]
    t_M = M = np.zeros((n, m), dtype=np.float32)
    ej = np.zeros((m,), dtype=np.float32)

    slice_obj = [slice(None, None, -1)] * len(filtertaps.shape)
    sliced_filtertaps_view = filtertaps[slice_obj]

    ps = ej.shape[0] + sliced_filtertaps_view.shape[0] - 1
    t_in1zpadded = in1zpadded = np.zeros(ps, ej.dtype)
    out = np.empty(ps, ej.dtype)

    for rowidx in range(n):
        for colidx in range(m):
            t_in1zpadded[colidx] = 1

            t_convolution = _correlateND(in1zpadded, sliced_filtertaps_view, out, 2)

            for kidx in range(k):
                t_M[rowidx, colidx] += t_convolution[<int>t_keep_indices[rowidx, kidx]] * t_V[rowidx, kidx]

            t_in1zpadded[colidx] = 0

    return M

它的速度快了 10 倍以上，但没有您的天方夜谭估计的那么快。话又说回来，这个计算一开始有点假；）。