NumExpr：加速Numpy、Pandas数学运算新利器！

量化投资与机器学习微信公众号

发布于 2020-06-29 16:07:22

2.5K0

发布于 2020-06-29 16:07:22

文章被收录于专栏：量化投资与机器学习量化投资与机器学习

作者：Sarkar 编译：1+1=6

前言

Numpy 和 Pandas 可能是用于数据科学（DS）和机器学习（ML）任务的两个最广泛使用的核心Python库。毋庸置疑，计算数值表达式的速度对于这些DS/ML任务至关重要，这两个库在这方面不会令人失望。

今天，我们又要给大家推荐一款利器：NumExpr。用来提高由Numpy和Pandas所产生的数学运算速度。

NumExpr的运算机制是怎么样的呢？

通常，表达式是使用 Python 编译函数编译的，提取变量并构建解析树结构。然后，这个树被编译成一个字节码程序，该程序使用所谓的“向量寄存器”（每个4096个元素宽）来描述基于元素的操作流。提高速度的关键是Numexpr一次处理元素块的能力。

它跳过了Numpy使用临时数组的做法，因为临时数组会浪费内存，而且对于大型数组，甚至无法装入缓存内存中。

另外，虚拟机完全是用C编写的，这使得它比本机Python更快。它也是多线程的，允许在合适的硬件上更快地并行化操作。

NumExpr支持在表达式中使用大量的数学运算符，但不支持条件运算符，如 if 或 else。

你也可以通过设置环境变量 NUMEXPR_MAX_THREAD 来控制你想要生成的线程的数量，以便用大型数组进行并行操作。目前，最大可能的线程数是64个，但是如果线程数高于底层CPU节点上可用的虚拟核数，就没有什么实际好处了。

下面是官方给出的解释：

GitHub地址：https://github.com/pydata/numexpr。

安装NumExpr库

和其他库一样：

pip install numexpr

根据源代码，“NumExpr是NumPy的快速数值表达式求值器。使用它，对数组进行操作的表达式可以得到加速，并且比在Python中进行相同的计算使用更少的内存。此外，它的多线程功能可以使用所有的内核——这通常会导致与NumPy相比性能的大幅提升。”(来源)

下面是具体文档，大家可以自行查看：

https://numexpr.readthedocs.io/projects/NumExpr3/en/latest/

标量-向量运算

从简单的数学运算开始。向Numpy数组添加一个标量，比如1。为了使用NumExpr包，我们所要做的就是将相同的计算包装在符号表达式中的特殊方法evaluate下：

a = np.arange(1e6)
b = np.arange(1e6)

%%timeit -n200 -r10
c = a+1
3.55 ms ± 52.1 µs per loop (mean ± std. dev. of 10 runs, 200 loops each)

%%timeit -n200 -r10
c = ne.evaluate("a + 1")
1.94 ms ± 86.5 µs per loop (mean ± std. dev. of 10 runs, 200 loops each)

速度有显著提升：从3.55ms提高到1.94ms！

两个数组运算

如下：

%%timeit -n100 -r10
c = 2*a+3*b
11.7 ms ± 177 µs per loop (mean ± std. dev. of 10 runs, 100 loops each)

%%timeit -n100 -r10
c = ne.evaluate("2*a+3*b")
2.14 ms ± 130 µs per loop (mean ± std. dev. of 10 runs, 100 loops each)

平均而言，计算时间从11.7 ms提高到了2.14 ms。

多数组复杂运算

让我们更进一步，在一个复杂的有理函数表达式中加入更多的数组。假设，我们想计算下面涉及5个Numpy数组的值，每个数组都有100万个随机数（从正态分布抽取）：

c = \frac{{a_1^2+2a_2+(3/a_3)}}{{\sqrt{a_4^2+a_5^2}}}

我们创建一个形状（1000000，5）的Numpy数组，并从中提取5个向量（1000000，1）用于有理函数：

a = np.random.normal(size=(1000000,5))
a1,a2,a3,a4,a5 = a[:,0],a[:,1],a[:,2],a[:,3],a[:,4]

%%timeit -n100 -r10
c = (a1**2+2*a2+(3/a3))/(np.sqrt(a4**2+a5**2))
47 ms ± 220 µs per loop (mean ± std. dev. of 10 runs, 100 loops each)

%%timeit -n100 -r10
ne.evaluate("(a1**2+2*a2+(3/a3))/(sqrt(a4**2+a5**2))")
3.96 ms ± 218 µs per loop (mean ± std. dev. of 10 runs, 100 loops each)

巨大的速度提升！直接从47ms下降到4ms。实际上，这是一个趋势，你会观察到：表达式变得越复杂，涉及的数组越多，使用Numexpr的速度提升就越快！

逻辑表达式 / bool过滤

我们并不局限于简单的算术表达式。Numpy数组最有用的特征之一是直接在包含逻辑运算符（如>或<）的表达式中使用它们来创建布尔过滤器或掩码。

我们可以用NumExpr做同样的操作，并加快过滤过程。我们检测欧氏距离测量涉及的4个向量是否大于某个阈值：

x1 = np.random.random(1000000)
x2 = np.random.random(1000000)
y1 = np.random.random(1000000)
y2 = np.random.random(1000000)

%%timeit -n100 -r10
c = np.sqrt((x1-x2)**2+(y1-y2)**2) > 0.5
23.2 ms ± 143 µs per loop (mean ± std. dev. of 10 runs, 100 loops each)

%%timeit -n100 -r10
c = ne.evaluate("sqrt((x1-x2)**2+(y1-y2)**2) > 0.5")
1.86 ms ± 112 µs per loop (mean ± std. dev. of 10 runs, 100 loops each)

%%timeit -n100 -r10
c = ne.evaluate("2*a+3*b > 3.5",optimization='moderate')
763 µs ± 85.4 µs per loop (mean ± std. dev. of 10 runs, 100 loops each)

在数据科学、机器学习pipeline中，这种过滤操作经常出现，你可以使用NumExpr表达式有策略地替换Numpy计算，这样可以节省很多计算时间。

复数

NumExpor也可以很好地处理复数，Python和Numpy本身就支持复数。这里有一个例子：

a = np.random.random(1000000)
b = np.random.random(1000000)

cplx = a + b*1j

%%timeit -n100 -r10
c = np.log10(cplx)
55.9 ms ± 159 µs per loop (mean ± std. dev. of 10 runs, 100 loops each)

%%timeit -n100 -r10
c = ne.evaluate("log10(cplx)")
9.9 ms ± 117 µs per loop (mean ± std. dev. of 10 runs, 100 loops each)

数组大小的影响

接下来，我们研究Numpy数组的大小对速度改进的影响。为此，我们选择一个简单的条件表达式，其中包含2*a+3*b < 3.5这样的两个数组，并绘制各种大小的相对执行时间（平均运行10次之后）。结果如下：

from time import time

result_np = {'Size':[],'Time':[]}
for i in [int(10**(j/5)) for j in range(25,40)]:
    a = np.random.random(size=i)
    b = np.random.random(size=i)
    times = [0]*10
    for j in range(10):
        t1 = time()
        c = (2*a+3*b > 3.5)
        t2 = time()
        times[j]=(t2-t1)*1000
    times = np.array(times)
    result_np['Size'].append(i)
    result_np['Time'].append(times.mean())

result_ne = {'Size':[],'Time':[]}
for i in [int(10**(j/5)) for j in range(25,40)]:
    a = np.random.random(size=i)
    b = np.random.random(size=i)
    times = [0]*10
    for j in range(10):
        t1 = time()
        c = ne.evaluate("2*a+3*b > 3.5")
        t2 = time()
        times[j]=(t2-t1)*1000
    times = np.array(times)
    result_ne['Size'].append(i)
    result_ne['Time'].append(times.mean())

def speed_benchmark(result1,result2,leg_text):
    """
    Plots timing results
    """
    plt.semilogx(result1['Size'],result1['Time'],c='blue',marker='o')
    plt.semilogx(result2['Size'],result2['Time'],c='k',marker='^')
    plt.grid(True)
    plt.legend(leg_text,fontsize=14)
    plt.xticks(fontsize=13)
    plt.yticks(fontsize=13)
    plt.xlabel('Number of elements in the array',fontsize=15)
    plt.ylabel("Time (milliseconds)",fontsize=15)
    plt.show()

speed_benchmark(result_np,result_ne,leg_text=['Just NumPy','With numexpr'])

pandas eval方法

这是一个对Python符号表达式（作为字符串）求值的Pandas方法。默认情况下，它使用NumExpr引擎来实现显著的加速：

https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.eval.html

使用以下代码，我们做一个简单的例子：构造四个DataFrame，每个数据包含50000行和100列（均匀随机数），并计算了一个涉及这些 DataFrames 的非线性变换。在一种情况下使用Pandas表达式，在另一种情况下使用pd.eval()方法。

nrows, ncols = 50000, 100

df1, df2, df3, df4 = [pd.DataFrame(np.random.randn(nrows, ncols)) for _ in range(4)]

%%timeit -n20 -r10
c=2*df1 - (df2/2) + (df3/df4)
55.8 ms ± 1.8 ms per loop (mean ± std. dev. of 10 runs, 20 loops each)

%%timeit -n20 -r10
pd.eval('2*df1 - (df2/2) + (df3/df4)')
17.3 ms ± 539 µs per loop (mean ± std. dev. of 10 runs, 20 loops each)

DataFrame大小的影响

我们对DataFrame的大小（行数，同时保持列数固定：100）对速度改进的影响进行了类似的分析。结果如下：

cols=100
result_no_eval = {'Size':[],'Time':[]}
for i in [int(10**(j/5)) for j in range(15,32)]:
    df1, df2, df3, df4 = [pd.DataFrame(np.random.randn(i, ncols)) for _ in range(4)]
    times = [0]*10
    for j in range(10):
        t1 = time()
        c = df1+df2+df3+df4
        t2 = time()
        times[j]=(t2-t1)*1000
    times = np.array(times)
    result_no_eval['Size'].append(i)
    result_no_eval['Time'].append(times.mean())

ncols=100
result_eval = {'Size':[],'Time':[]}
for i in [int(10**(j/5)) for j in range(15,32)]:
    df1, df2, df3, df4 = [pd.DataFrame(np.random.randn(i, ncols)) for _ in range(4)]
    times = [0]*10
    for j in range(10):
        t1 = time()
        c = ne.evaluate("df1+df2+df3+df4")
        t2 = time()
        times[j]=(t2-t1)*1000
    times = np.array(times)
    result_eval['Size'].append(i)
    result_eval['Time'].append(times.mean())

def speed_benchmark_pd(result1,result2,leg_text):
    """
    Plots timing results
    """
    plt.semilogx(result1['Size'],result1['Time'],c='blue',marker='o')
    plt.semilogx(result2['Size'],result2['Time'],c='k',marker='^')
    plt.grid(True)
    plt.legend(leg_text,fontsize=14)
    plt.xticks(fontsize=13)
    plt.yticks(fontsize=13)
    plt.xlabel('Number of rows in the DataFrame',fontsize=15)
    plt.ylabel("Time (milliseconds)",fontsize=15)
    plt.show()
    
speed_benchmark_pd(result_no_eval,result_eval,
                   leg_text=['Normal boring Pandas','With pd.eval'])