专栏首页量化投资与机器学习NumExpr:加速Numpy、Pandas数学运算新利器!

NumExpr:加速Numpy、Pandas数学运算新利器!

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

1

前言

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

2

安装NumExpr库

和其他库一样:

pip install numexpr

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

下面是具体文档,大家可以自行查看:

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

3

标量-向量运算

从简单的数学运算开始。向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!

4

两个数组运算

如下:

%%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

多数组复杂运算

让我们更进一步,在一个复杂的有理函数表达式中加入更多的数组。假设,我们想计算下面涉及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的速度提升就越快!

6

逻辑表达式 / 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计算,这样可以节省很多计算时间。

7

复数

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)

8

数组大小的影响

接下来,我们研究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'])

9

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)

10

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

NumExpr可能不适用与所有操作,但是大部分的数据科学、统计建模等都有应用之处。而且对代码的修改也很小。希望大家可以有所收获!

本文分享自微信公众号 - 量化投资与机器学习(Lhtz_Jqxx),作者:QIML编辑部

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-06-23

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 最新 | 基于回声状态网络预测股票价格(附代码)

    “There (is) order and even great beauty in what looks like total chaos. If we lo...

    量化投资与机器学习微信公众号
  • 【ML系列】一招鲜,判断哪些输入特征对神经网络是重要的!

    没错,但这不是今天的重点。我们想知道的是输入特征对神经网络的预测计算有多重要。例如,通过学习时间、年龄、身高和缺席人数等几个预测因素来预测谁会通过考试。直觉上,...

    量化投资与机器学习微信公众号
  • 神经网络在算法交易上的应用系列——多元时间序列

    之前的文章已经介绍了几种预测时间序列的方法:如何规范化数据,以实值或二进制变量的形式进行预测,以及如何处理高噪声中的过拟合。在上一篇文章中,我们只用了经过一些转...

    量化投资与机器学习微信公众号
  • 衡量直播平台的推流效果,主要看这5个指标

    近年来,网络直播呈现爆发式增长,上百家平台,超百亿规模,3亿多用户,上市公司和明星企业崛起,俨然成为产业。光环加持的背后,腾讯云直播TLive平台,抓住机遇、持...

    腾讯云视频
  • 数据结构之链表

    1、线性数据结构,动态数组、栈、队列,底层依托静态数组,靠resize解决固定容量问题。

    别先生
  • Python教程01:词云Word Cloud微服务

    在日常项目中要生成类似上图这样的云词图片. 这里有一个python库可以完美的解决生成云词图片 https://github.com/amueller/word...

    mojocn
  • JAVA NIO Socket通道

    DatagramChannel和SocketChannel都实现定义读写功能,ServerSocketChannel不实现,只负责监听传入的连接,并建立新的So...

    WindWant
  • 【Spark篇】---Spark故障解决(troubleshooting)

    1) connection timeout ----shuffle file cannot find

    LhWorld哥陪你聊算法
  • 手把手教你使用PCA进行数据降维

    对数据降维可以帮助我们提取数据集的主要信息,即将原始的高维特征空间压缩到低纬度的特征子空间。数据降维是用于提高计算效率的典型手段,另一个好处是也能够减小维度诅咒...

    HuangWeiAI
  • 小程序formid埋点

    只有form提交的时候才会有formId,这样用户的formId可能数量比较少,不能实现发送很多微信模板消息。所以要用伪装表单的方式来实现获取formId。

    薛定喵君

扫码关注云+社区

领取腾讯云代金券