NumPy比Numba和Cython更快，如何改进Numba代码

Question

我这里有一个简单的例子，帮助我了解如何使用numba和cython。我对numba和cython都很新。我已经尝试了最好的方法来使numba快速，并且在某种程度上，对于cython也是如此。但是对于float64，我的numpy代码几乎比numba快2倍，如果使用float32，则快2倍以上。不确定我错过了什么。

我想也许问题不再是编码方面，而更多地涉及编译器等，这些我不太熟悉。

我已经浏览了许多有关numpy、numba和cython的stackoverflow帖子，没有直接的答案。

numpy版本：

def py_expsum(x):
    return np.sum( np.exp(x) )

Numba版本：

@numba.jit( nopython=True)    
def nb_expsum(x):
    nx, ny = x.shape
    val = 0.0
    for ix in range(nx):
        for iy in range(ny):
            val += np.exp(x[ix, iy])
    return val

Cython版本：

import numpy as np
import cython
from libc.math cimport exp

@cython.boundscheck(False) 
@cython.wraparound(False)
cpdef double cy_expsum2 ( double[:,:] x, int nx, int ny ):
    cdef: 
        double val = 0.0
        int ix, iy    
    for ix in range(nx):
        for iy in range(ny):
            val += exp(x[ix, iy])
    return val

玩弄大小为2000 x 1000的数组，并循环100次。对于numba，第一次启用不计入循环。

使用Python 3（Anaconda发行版），Windows 10。

               float64       /   float32
    1. numpy : 0.56 sec      /   0.23 sec
    2. numba : 0.93 sec      /   0.74 sec      
    3. cython: 0.83 sec

cython几乎与numba相似。但对我来说，一个大问题是为什么numba跑不过numpy? 我做错了什么或者遗漏了什么？其他因素如何影响性能，我该如何找出原因？

- Ong Beng Seong

---

你应该使用 math.exp 而不是 np.exp。 - Divakar

什么是打字错误？什么是相同的信息？ - Divakar

错别字已更正。math.exp 没有帮助。 - Ong Beng Seong

Numpy 可能会并行处理指数运算。虽然你也可以使用 Cython（和可能的 Numba）来实现这一点，但你很可能无法显著超过 Numpy。为什么不直接使用 Numpy 呢？ - DavidW

很难超越numpy向量化代码的性能。但是如果你想要一点性能提升，可以使用numexpr，例如：ne.evaluate('sum(exp(x))')。 - Brenlla

Answer 1

如我们所见，行为取决于使用哪个numpy分布。

本答案将重点介绍使用Intel的VML（向量数学库）的Anaconda分布，但在其他硬件和numpy版本下，效果可能有所不同。

还将展示如何通过Cython或numexpr利用VML，在不使用Anaconda分布的情况下进行一些numpy运算，因为Anaconda分布在底层为某些numpy操作插入了VML。

---

我可以在以下维度上复现您的结果。

N,M=2*10**4, 10**3
a=np.random.rand(N, M)

我得到：

%timeit py_expsum(a)  #   87ms
%timeit nb_expsum(a)  #  672ms
%timeit nb_expsum2(a)  #  412ms

计算时间的绝大部分（约90%）用于exp函数的评估，正如我们将看到的，这是一个CPU密集型任务。

快速查看top统计数据显示，numpy的版本已经执行了并行化，但numba没有。然而，在我的只有两个处理器的虚拟机上，仅靠并行化就无法解释巨大的7倍差异（如DavidW的nb_expsum2版本所示）。

通过perf为两个版本编写代码进行分析，结果如下：

nb_expsum

Overhead  Command  Shared Object                                      Symbol                                                             
  62,56%  python   libm-2.23.so                                       [.] __ieee754_exp_avx
  16,16%  python   libm-2.23.so                                       [.] __GI___exp
   5,25%  python   perf-28936.map                                     [.] 0x00007f1658d53213
   2,21%  python   mtrand.cpython-37m-x86_64-linux-gnu.so             [.] rk_random

py_expsum

  31,84%  python   libmkl_vml_avx.so                                  [.] mkl_vml_kernel_dExp_E9HAynn                                   ▒
   9,47%  python   libiomp5.so                                        [.] _INTERNAL_25_______src_kmp_barrier_cpp_38a91946::__kmp_wait_te▒
   6,21%  python   [unknown]                                          [k] 0xffffffff8140290c                                            ▒
   5,27%  python   mtrand.cpython-37m-x86_64-linux-gnu.so             [.] rk_random  

可以看到：numpy在其内部使用了Intel的并行向量化mkl/vml版本，它比numba使用的gnu-math-library (lm.so)，或者是numba的并行版本或者是cython版本都更快。通过使用并行化可以稍微平衡一下，但仍然是mkl的向量化版本优于numba和cython。

然而，仅仅看一个大小的性能并不能很好地说明问题，对于exp(和其他超越函数)，需要考虑两个因素：

• 数组中元素的数量——缓存效应和不同大小的不同算法（在numpy中并不罕见）可能导致不同的性能。
• 根据x-值，需要计算exp(x)的时间也不同。通常有三种不同类型的输入会导致不同的计算时间：非常小、正常和非常大（具有非有限结果）

我使用perfplot来可视化结果（请参见附录中的代码）。对于“正常”范围，我们获得以下性能：

虽然0.0的性能类似，但我们可以看到，一旦结果变为无穷大，Intel的VML会受到相当大的负面影响：

然而还有其他要观察到的事情：

• 对于向量大小<=8192 = 2^13，numpy使用未并行化的glibc版本的exp（numba和cython也使用相同的版本）。
• 我使用的Anaconda分发版覆盖了numpy的功能并插入了Intel的VML库，用于大小>8192的向量化和并行化 - 这解释了在大小约为10^4的情况下运行时间下降。
• numba轻松击败了通常的glibc版本（numpy开销太大），但对于更大的数组来说，如果numpy不切换到VML，差别不会很大。
• 它似乎是一个CPU绑定的任务——我们无法在任何地方看到缓存边界。
• 并行化的numba版本只有在有500个以上的元素时才有意义。

那么，有什么后果呢？

1. 如果不超过8192个元素，应该使用numba版本。
2. 否则，应该使用numpy版本（即使没有VML插件也不会失去太多）。

NB：numba不能自动使用Intel的VML中的vdExp（在评论中部分建议），因为它单独计算exp(x)，而VML操作整个数组。

---

当写入和加载数据时，可以通过numpy版本使用以下算法来减少缓存未命中：

1. 对适合缓存但也不太小（开销较大）的一部分数据执行VML的vdExp。
2. 对

   %%cython -L=<path_mkl_libs> --link-args=-Wl,-rpath=<path_mkl_libs> --link-args=-Wl,--no-as-needed -l=mkl_intel_ilp64 -l=mkl_core -l=mkl_gnu_thread -l=iomp5
   # path to mkl can be found via np.show_config()
   # which libraries needed: https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor
   
   # another option would be to wrap mkl.h:
   cdef extern from *:
       """
       // MKL_INT is 64bit integer for mkl-ilp64
       // see https://software.intel.com/en-us/mkl-developer-reference-c-c-datatypes-specific-to-intel-mkl
       #define MKL_INT long long int
       void  vdExp(MKL_INT n, const double *x, double *y);
       """
       void vdExp(long long int n, const double *x, double *y)
   
   def cy_expsum(const double[:,:] v):
           cdef:
               double[1024] w;
               int n = v.size
               int current = 0;
               double res = 0.0
               int size = 0
               int i = 0
           while current<n:
               size = n-current
               if size>1024:
                   size = 1024
               vdExp(size, &v[0,0]+current, w)
               for i in range(size):
                   res+=w[i]
               current+=size
           return res
   

   然而，这正是numexpr所做的，它也使用Intel的vml作为后端:

    import numexpr as ne
    def ne_expsum(x):
        return ne.evaluate("sum(exp(x))")
   

   关于时间，我们可以看到如下：

   

   以下是值得注意的细节:

   • 对于较大的数组，numpy、numexpr 和 cython 版本几乎具有相同的性能 - 这并不令人惊讶，因为它们使用相同的 vml 功能。
   • 在这三个选项中，cython 版本的开销最小，而 numexpr 版本的开销最大。
   • 由于不是每个 numpy 发布版本都包含 mvl 功能，因此 numexpr 版本可能是最容易编写的。

   ---

   代码清单:

   图表:

   import numpy as np
   def py_expsum(x):
       return np.sum(np.exp(x))
   
   import numba as nb
   @nb.jit( nopython=True)    
   def nb_expsum(x):
       nx, ny = x.shape
       val = 0.0
       for ix in range(nx):
           for iy in range(ny):
               val += np.exp( x[ix, iy] )
       return val
   
   @nb.jit( nopython=True, parallel=True)    
   def nb_expsum2(x):
       nx, ny = x.shape
       val = 0.0
       for ix in range(nx):
           for iy in nb.prange(ny):
               val += np.exp( x[ix, iy]   )
       return val
   
   import perfplot
   factor = 1.0 # 0.0 or 1e4
   perfplot.show(
       setup=lambda n: factor*np.random.rand(1,n),
       n_range=[2**k for k in range(0,27)],
       kernels=[
           py_expsum, 
           nb_expsum,
           nb_expsum2, 
           ],
       logx=True,
       logy=True,
       xlabel='len(x)'
       )
   

- ead

---

非常感谢您的建议。我不知道numpy正在进行并行化处理。因此，更公平的测试将是强制numba和cython也进行并行化处理。 - Ong Beng Seong

你的结果看起来像是numpy使用了Intel SVML，而numba和cython没有使用。SVML可以很容易地安装。https://numba.pydata.org/numba-doc/dev/user/performance-tips.html - max9111

@max9111 我认为numba也无法使用MVL中的vdExp函数，因为它是对数组进行操作而不是单个值。 - ead

如果数组的值超过8192个，numpy的结果也可能会略微改变。https://dev59.com/-VMI5IYBdhLWcg3wyunE 你尝试安装SVML了吗？（conda install -c numba icc_rt） - max9111

我怎样在Linux和Windows中找到共享对象列表？我猜你说的是通过一个叫做perf的工具在Linux上实现的吧？ - Ong Beng Seong

Answer 2

添加并行处理。在Numba中，这只涉及将外部循环更改为prange并向jit选项中添加parallel=True：

@numba.jit( nopython=True,parallel=True)    
def nb_expsum2(x):
    nx, ny = x.shape
    val = 0.0
    for ix in numba.prange(nx):
        for iy in range(ny):
            val += np.exp( x[ix, iy]   )
    return val

在我的电脑上，这比非并行版本快了3.2倍。话虽如此，在我的电脑上，Numba和Cython都比Numpy更好。您还可以在Cython中进行并行处理 - 我没有在这里测试过，但我希望其性能与Numba类似。（还要注意的是，对于Cython，您可以从x.shape [0]和x.shape [1]获取nx和ny，因此您不必关闭边界检查，然后完全依赖用户输入来保持在范围内）。

- DavidW

---

谢谢DavidW。我不知道NumPy自动使用并行处理。我刚刚尝试了Numba的并行选项，但对我来说没有区别。至于你的情况，Numba/Cython击败了Numba。我发现在我的较慢的笔记本电脑上（也是较少的核心），有一个效果。这是击败NumPy的唯一原因吗？由于核心较少，NumPy的并行处理能力也较弱吗？ - Ong Beng Seong

确保您对Numba并行化的代码进行两个更改。很难知道相对速度的确切原因 - 它可能取决于编译器、CPU以及编译选项。然而，大体上有两件主要事情可以变化：是否并行运行以及创建临时数组（NumPy版本会创建临时数组，而其他版本则不会）。 - DavidW

Answer 3

这取决于经验实现和并行化

如果您在Numpy中使用Intel SVML，请在其他软件包（如Numba，Numexpr或Cython）中也使用它。 Numba性能提示

如果Numpy命令已经并行化，请尝试在Numba或Cython中并行化它。

代码

import os
#Have to be before importing numpy
#Test with 1 Thread against a single thread Numba/Cython Version and
#at least with number of physical cores against parallel versions
os.environ["MKL_NUM_THREADS"] = "1" 

import numpy as np

#from version 0.43 until 0.47 this has to be set before importing numba
#Bug: https://github.com/numba/numba/issues/4689
from llvmlite import binding
binding.set_option('SVML', '-vector-library=SVML')
import numba as nb

def py_expsum(x):
    return np.sum( np.exp(x) )

@nb.njit(parallel=False,fastmath=True) #set it to True for a parallel version  
def nb_expsum(x):
    val = nb.float32(0.)#change this to float64 on the float64 version
    for ix in nb.prange(x.shape[0]):
        for iy in range(x.shape[1]):
            val += np.exp(x[ix,iy])
    return val

N,M=2000, 1000
#a=np.random.rand(N*M).reshape((N,M)).astype(np.float32)
a=np.random.rand(N*M).reshape((N,M))

基准测试

#float64
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "1" 
#7.44 ms ± 86.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "6" 
#4.83 ms ± 139 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=false
#2.49 ms ± 25.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) ##parallel=true
#568 µs ± 45.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

#float32
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "1" 
#3.44 ms ± 66.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit py_expsum(a) #os.environ["MKL_NUM_THREADS"] = "6" 
#2.59 ms ± 35.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nb_expsum(a) #parallel=false
#1 ms ± 12.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit nb_expsum(a) #parallel=true
#252 µs ± 19.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

使用SVML进行性能绘图

import numpy as np

#from version 0.43 until 0.47 this has to be set before importing numba
#Bug: https://github.com/numba/numba/issues/4689
from llvmlite import binding
binding.set_option('SVML', '-vector-library=SVML')
import numba as nb

def py_expsum(x):
    return np.sum(np.exp(x))

@nb.jit( nopython=True,parallel=False,fastmath=False)    
def nb_expsum_single_thread(x):
    nx, ny = x.shape
    val = 0.0
    for ix in range(nx):
        for iy in range(ny):
            val += np.exp( x[ix, iy] )
    return val

#fastmath makes SIMD-vectorization possible 
#val+=some_value is not vectorizable (scalar depends on scalar)
#This would also prevents the usage of SVML
@nb.jit( nopython=True,parallel=False,fastmath=True)    
def nb_expsum_single_thread_vec(x):
    nx, ny = x.shape
    val = 0.0
    for ix in range(nx):
        for iy in range(ny):
            val += np.exp( x[ix, iy] )
    return val

@nb.jit(nopython=True,parallel=True,fastmath=False)    
def nb_expsum_parallel(x):
    nx, ny = x.shape
    val = 0.0
    #parallelization over the outer loop is almost every time faster
    #except for rare cases like this (x.shape -> (1,n))
    for ix in range(nx):
        for iy in nb.prange(ny):
            val += np.exp( x[ix, iy] )
    return val

#fastmath makes SIMD-vectorization possible 
#val+=some_value is not vectorizable (scalar depends on scalar)
#This would also prevents the usage of SVML
@nb.jit(nopython=True,parallel=True,fastmath=True)    
def nb_expsum_parallel_vec(x):
    nx, ny = x.shape
    val = 0.0
    #parallelization over the outer loop is almost every time faster
    #except for rare cases like this (x.shape -> (1,n))
    for ix in range(nx):
        for iy in nb.prange(ny):
            val += np.exp( x[ix, iy] )
    return val

import perfplot
factor = 1.0 # 0.0 or 1e4
perfplot.show(
    setup=lambda n: factor*np.random.rand(1,n),
    n_range=[2**k for k in range(0,27)],
    kernels=[
        py_expsum,
        nb_expsum_single_thread,
        nb_expsum_single_thread_vec,
        nb_expsum_parallel,
        nb_expsum_parallel_vec,
        cy_expsum
        ],
    logx=True,
    logy=True,
    xlabel='len(x)'
    )

检查是否使用了SVML

可以用来检查一切是否正常运行。

def check_SVML(func):
    if 'intel_svmlcc' in func.inspect_llvm(func.signatures[0]):
        print("found")
    else:
        print("not found")

check_SVML(nb_expsum_parallel_vec)
#found

- max9111