为什么 Cython 比矢量化 NumPy 慢？

Question

考虑以下 Cython 代码：

cimport cython
cimport numpy as np
import numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
def test_memoryview(double[:] a, double[:] b):
    cdef int i
    for i in range(a.shape[0]):
        a[i] += b[i]

@cython.boundscheck(False)
@cython.wraparound(False)
def test_numpy(np.ndarray[double, ndim=1] a, np.ndarray[double, ndim=1] b):
    cdef int i
    for i in range(a.shape[0]):
        a[i] += b[i]

def test_numpyvec(a, b):
    a += b

def gendata(nb=40000000):
    a = np.random.random(nb)
    b = np.random.random(nb)
    return a, b

在解释器中运行它会产生（经过几次运行以预热缓存）：

In [14]: %timeit -n 100 test_memoryview(a, b)
100 loops, best of 3: 148 ms per loop

In [15]: %timeit -n 100 test_numpy(a, b)
100 loops, best of 3: 159 ms per loop

In [16]: %timeit -n 100 test_numpyvec(a, b)
100 loops, best of 3: 124 ms per loop

# See answer below :
In [17]: %timeit -n 100 test_raw_pointers(a, b)
100 loops, best of 3: 129 ms per loop

我尝试了不同的数据集大小，并且始终让矢量化 NumPy 函数比编译后的 Cython 代码运行得更快，而我期望 Cython 在性能方面与矢量化 NumPy 相当。

我是否忘记了 Cython 代码中的优化？ NumPy 是否使用某些东西（BLAS？）来使这些简单的操作运行得更快？我可以提高这段代码的性能吗？

更新：原始指针版本似乎与 NumPy 相当。因此，使用内存视图或 NumPy 索引显然存在一些开销。

Answer 1

在我的机器上差异不是很大，但我几乎可以通过像这样更改 numpy 和内存视图函数来消除它

@cython.boundscheck(False)
@cython.wraparound(False)
def test_memoryview(double[:] a, double[:] b):
    cdef int i, n=a.shape[0]
    for i in range(n):
        a[i] += b[i]

@cython.boundscheck(False)
@cython.wraparound(False)
def test_numpy(np.ndarray[double] a, np.ndarray[double] b):
    cdef int i, n=a.shape[0]
    for i in range(n):
        a[i] += b[i]

然后，当我从 Cython 编译 C 输出时，我使用标志 -O3 和 -march=native。 这似乎表明时序的差异来自于使用不同的编译器优化。

我使用 64 位版本的 MinGW 和 NumPy 1.8.1。 您的结果可能会因您的软件包版本、硬件、平台和编译器而异。

如果您使用的是 IPython notebook 的 Cython 魔法，您可以通过将 %%cython 替换为 %%cython -f -c=-O3 -c=-march=native 来强制使用其他编译器标志进行更新

如果您为 cython 模块使用标准 setup.py，则可以在创建传递给 distutils.setup 的扩展对象时指定 extra_compile_args 参数。

注意：我在指定 NumPy 数组的类型时删除了 ndim=1 标志，因为它不是必需的。 无论如何，该值默认为 1。

Answer 2

稍微提高速度的一个变化是指定步幅：

def test_memoryview_inorder(double[::1] a, double[::1] b):
    cdef int i
    for i in range(a.shape[0]):
        a[i] += b[i]

Answer 3

另一种选择是使用原始指针（以及避免重复@cython... 的全局指令）：

#cython: wraparound=False
#cython: boundscheck=False
#cython: nonecheck=False

#...

cdef ctest_raw_pointers(int n, double *a, double *b):
    cdef int i
    for i in range(n):
        a[i] += b[i]

def test_raw_pointers(np.ndarray[double, ndim=1] a, np.ndarray[double, ndim=1] b):
    ctest_raw_pointers(a.shape[0], &a[0], &b[0])