numpy比numba和cython快，如何改进numba代码

Question

这里我有一个简单的例子来帮助我理解使用numba和cython。我是numba和cython的新手。我已经尽力将所有的技巧结合起来，使numba快速，在某种程度上，对于cython也是如此，但是我的numpy代码几乎比numba ( float64)快2倍，如果使用float32，则比numba快2倍。不知道我在这里错过了什么。

我在想，也许问题不再是编码，而是更多关于编译器之类的问题，这一点我不太熟悉。

我浏览了很多关于numpy，numba和cython的堆叠文章，却没有找到任何直接的答案。

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分发版)，窗口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的运行时？我在这里做错了什么或者错过了什么？其他因素如何促成，我如何找出呢？

Answer 1

正如我们将要看到的，行为取决于使用numpy-分发的是什么。

这个答案将集中在英特尔的VML (矢量数学库)的Anacoda发行版上，在另一个硬件和numpy版本的情况下，磨耗可能会有所不同。

它还将展示如何通过Cython或numexpr来使用VML，如果不使用Anacoda-分发，它将在VML中插入一些numpy操作。

我可以复制您的结果，如下所示

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-statistics显示，numpy的版本已被平分化，但numba的情况并非如此。然而，在只有两个处理器的VM上，并行化本身并不能解释因子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在引擎盖下使用英特尔的Paralliized矢量化mkl/vml版本，它的性能轻松超过numba使用的gnu-数学库(lm.so)版本(或者是numba的并行版本，或者是cython的版本)。我们可以通过使用分段技术来平平地面，但是mkl的矢量化版本的性能仍然优于numba和cython。

然而，只看到一个大小的性能并不是很有启发性，对于exp (对于其他先验函数)，有两个维度需要考虑：

• 数组缓存效应中的元素数和不同大小的不同算法(在numpy中不是闻所未闻)可以导致不同的性能。
• 根据x-value，计算exp(x)需要不同的时间。通常，有三种不同类型的输入，导致不同的计算时间:非常小，正常和非常大(与非有限的结果)。

我正在使用完善图来可视化结果(参见附录中的代码)。对于“正常”范围，我们得到以下性能：

​

​

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

​

​

然而，还有其他需要注意的事情：

• 对于向量大小，<= 8192 = 2^13 numpy使用非并行化的glibc版本的exp (同样的numba和cython也在使用)。
• Anaconda发行版，我使用的是> 8192的覆盖numpy的功能并插入Intel的VML库，它是矢量化和并行化的--这解释了大小10^4的运行时间减少的原因。
• numba轻松地击败了通常的glibc版本(对于numpy来说开销太大)，但是对于更大的数组(如果numpy不会切换到VML)则不会有太大的差别。
• 它似乎是一个CPU绑定的任务--我们看不到缓存--任何地方的边界。
• 只有当有超过500个元素时，Parallized版本才有意义。

那么后果是什么呢？

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

注: numba不能自动使用英特尔VML中的vdExp (如注释中部分建议的那样)，因为它单独计算exp(x)，而VML则在整个数组上运行。

在写入和加载数据时，可以减少缓存丢失，该数据由numpy版本使用以下算法执行：

1. 对适合高速缓存但也不太小(开销)的部分数据执行VML的vdExp。
2. 总结得到的工作数组。
3. 为下一部分数据执行1.+2。直到整个数据被处理为止。

然而，与numpy的版本相比，我不会期望获得超过10% (但也许我错了)，因为无论如何，90%的计算时间都是用在MVL上的。

然而，下面是Cython中一个可能的快速和肮脏的实现：

%%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所做的，它也使用英特尔的vml作为后端：

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

至于时间安排，我们可以看到以下情况：

​

​

以下是值得注意的细节：

• numpy、numexpr和cython版本对于更大的数组具有几乎相同的性能--这并不奇怪，因为它们使用相同的vml-功能。
• 从这三个版本中，cython版本的开销最小，numexpr的开销最大。
• numexpr-版本可能是最容易编写的(考虑到并非每个numpy发行版plugsin mvl-功能)。

清单：

情节：

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)'
    )

Answer 2

添加并行化。在Numba中，这只需要将外部循环prange和parallel=True添加到jit选项中：

@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

在我的PC上，它比非并行版本加速了3.2倍。这就是说，在我的个人电脑上，Numba和Cython都像写的那样击败了Numpy。

您也可以做Cython中的并行化 --我还没有在这里测试过它，但是我希望在性能上与Numba类似。(还请注意，对于Cython，您可以从nx和x.shape[1]获得x.shape[0]和x.shape[1]，这样您就不必关闭边界检查，然后完全依赖用户输入才能保持在界限内)。

Answer 3

它依赖于exp的实现和并行化。

如果在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