Numpy为什么在这个cython例程中胜过3倍

Question

我刚开始尝试使用cython,作为第一个练习,我创建了一个函数的以下(重新)实现来计算数组中每个元素的sin.所以这是我的sin.pyx

from numpy cimport ndarray, float64_t
import numpy as np

cdef extern from "math.h":
    double sin(double x)

def sin_array(ndarray[float64_t, ndim=1] arr):
    cdef int n = len(arr)
    cdef ndarray h = np.zeros(n, dtype=np.float64)
    for i in range(n):
        h[i] = sin(arr[i])
    return h

我还为此创建了以下setup.py

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

import numpy

ext = Extension("sin", sources=["sin.pyx"])

setup(ext_modules=[ext],
      cmdclass={"build_ext": build_ext},
      include_dirs=[numpy.get_include()])

所以这会创建我的*.so文件.我将其导入python并创建1000个随机数,例如

import sin
import numpy as np

x = np.random.randn(1000)

%timeit sin.sin_array(x)
%timeit np.sin(x)

Numpy赢了3倍.为什么？我认为对输入数组的类型和维度做出非常明确的假设的函数在这里可以更具竞争力.当然,我也明白numpy是不可思议的聪明,但很有可能我在这里做了一些愚蠢的事......

请注意,本练习的目的是不重写更快的sin函数,而是为我们的一些内部工具创建一些cython包装器,但这是以后的另一个问题......

Answer 1

Cython的注释功能,cython -a filename.pyx是你的朋友.它会生成一个html文件,您可以在浏览器中加载它,并突出显示未经过优化的代码行.您可以单击一行以查看生成的c代码.

在这种情况下,问题似乎h是没有正确输入.如果您只是键入一个数组,因为ndarray您告诉Cython它是一个数组,但是您没有给cython足够的信息来告诉它如何有效地索引它,您必须提供类型和形状信息.您已在函数声明中正确完成此操作.

我想,一旦修复后,性能将具有可比性,但如果没有注释会告诉你什么是错的.如果cython仍然较慢,那么numpy可能使用比标准c更快的sin函数(你可以获得更快的sin近似值,如果感兴趣的话可以使用谷歌搜索).

Answer 2

以下是使用ipython中的cython magic在我的机器上的几个变体和性能(可能会有所不同):

%%cython --compile-args=-O3 -a

import numpy as np
cimport numpy as np
import cython

from libc.math cimport sin

def sin_array(np.ndarray[np.float64_t, ndim=1] arr):
    cdef int n = len(arr)
    cdef np.ndarray h = np.zeros(n, dtype=np.float64)
    for i in range(n):
        h[i] = sin(arr[i])
    return h

@cython.boundscheck(False)
@cython.wraparound(False)
def sin_array1(np.ndarray[np.float64_t, ndim=1] arr):
    cdef int n = arr.shape[0]
    cdef unsigned int i
    cdef np.ndarray[np.float64_t, ndim=1] h = np.empty_like(arr)
    for i in range(n):
        h[i] = sin(arr[i])
    return h


@cython.boundscheck(False)
@cython.wraparound(False)
def sin_array2(np.float64_t[:] arr):
    cdef int n = arr.shape[0]
    cdef unsigned int i
    cdef np.ndarray[np.float64_t, ndim=1] h = np.empty(n, np.float64)
    cdef np.float64_t[::1] _h = h
    for i in range(n):
        _h[i] = sin(arr[i])
    return h

为了踢,我投入了一个Numba jitted方法:

import numpy as np
import numba as nb

@nb.jit
def sin_numba(x):
    n = x.shape[0]
    h = np.empty(n, np.float64)
    for k in range(n):
        h[k] = np.sin(x[k])

    return h

时间安排:

In [25]:

x = np.random.randn(1000)

%timeit np.sin(x)
%timeit sin_array(x)
%timeit sin_array1(x)
%timeit sin_array2(x)
%timeit sin_numba(x)
10000 loops, best of 3: 27 µs per loop
10000 loops, best of 3: 80.3 µs per loop
10000 loops, best of 3: 28.7 µs per loop
10000 loops, best of 3: 32.8 µs per loop
10000 loops, best of 3: 31.4 µs per loop

numpy内置仍然是最快的(但只是一点点),考虑到不指定任何类型信息的简单性,numba性能非常好.

更新:

看看各种阵列尺寸也总是好的.以下是10000个元素数组的时序:

In [26]:

x = np.random.randn(10000)

%timeit np.sin(x)
%timeit sin_array(x)
%timeit sin_array1(x)
%timeit sin_array2(x)
%timeit sin_numba(x)
1000 loops, best of 3: 267 µs per loop
1000 loops, best of 3: 783 µs per loop
1000 loops, best of 3: 267 µs per loop
1000 loops, best of 3: 268 µs per loop
1 loops, best of 3: 287 µs per loop

在这里,您可以看到原始方法和np.sin调用的优化版本之间几乎相同的时序,指向cython或返回中数据结构初始化的一些开销.在这些条件下,Numba的情况稍差.