为什么numpy的einsum比numpy的内置函数更快?

让我们从三个dtype=np.double数组开始。 计时是在英特尔CPU上使用nccu 1.7.1与icc编译并链接到intel的mkl 。 使用gcc而不使用mkl编译的numpy 1.6.1的AMD cpu也用于验证时序。 请注意,时序规模几乎与系统规模成线性关系,并不是因为numpy函数的小开销,因为if声明这些差异将以毫秒而非毫秒显示:

arr_1D=np.arange(500,dtype=np.double)
large_arr_1D=np.arange(100000,dtype=np.double)
arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500)
arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500)

首先让我们看一下np.sum函数:

np.all(np.sum(arr_3D)==np.einsum('ijk->',arr_3D))
True

%timeit np.sum(arr_3D)
10 loops, best of 3: 142 ms per loop

%timeit np.einsum('ijk->', arr_3D)
10 loops, best of 3: 70.2 ms per loop

鲍尔斯:

np.allclose(arr_3D*arr_3D*arr_3D,np.einsum('ijk,ijk,ijk->ijk',arr_3D,arr_3D,arr_3D))
True

%timeit arr_3D*arr_3D*arr_3D
1 loops, best of 3: 1.32 s per loop

%timeit np.einsum('ijk,ijk,ijk->ijk', arr_3D, arr_3D, arr_3D)
1 loops, best of 3: 694 ms per loop

外部产品:

np.all(np.outer(arr_1D,arr_1D)==np.einsum('i,k->ik',arr_1D,arr_1D))
True

%timeit np.outer(arr_1D, arr_1D)
1000 loops, best of 3: 411 us per loop

%timeit np.einsum('i,k->ik', arr_1D, arr_1D)
1000 loops, best of 3: 245 us per loop

以上所有都是np.einsum两倍。 这些应该是苹果比较的苹果,因为所有东西都是dtype=np.double 。 我希望在这样的操作中加快速度:

np.allclose(np.sum(arr_2D*arr_3D),np.einsum('ij,oij->',arr_2D,arr_3D))
True

%timeit np.sum(arr_2D*arr_3D)
1 loops, best of 3: 813 ms per loop

%timeit np.einsum('ij,oij->', arr_2D, arr_3D)
10 loops, best of 3: 85.1 ms per loop

无论axes选择如何, np.kronnp.innernp.outernp.kronnp.sum来说似乎至少快两倍。 主要的例外是np.dot因为它从BLAS库调用DGEMM。 那么为什么np.einsum比其他numpy函数更快呢?

DGEMM案件的完整性:

np.allclose(np.dot(arr_2D,arr_2D),np.einsum('ij,jk',arr_2D,arr_2D))
True

%timeit np.einsum('ij,jk',arr_2D,arr_2D)
10 loops, best of 3: 56.1 ms per loop

%timeit np.dot(arr_2D,arr_2D)
100 loops, best of 3: 5.17 ms per loop

主要理论来自@sebergs评论, np.einsum可以使用SSE2,但numpy的ufuncs不会在numpy 1.8之前(请参阅更改日志)。 我相信这是正确的答案,但一直未能证实。 通过改变输入数组的dtype并观察速度差异以及不是每个人都遵守相同的定时趋势这一事实,可以找到一些有限的证明。


首先,关于这个在numpy列表上有很多过去的讨论。 例如,请参阅:http://numpy-discussion.10968.n7.nabble.com/poor-performance-of-sum-with-sub-machine-word-integer-types-td41.html http:// numpy- discussion.10968.n7.nabble.com/odd-performance-of-sum-td3332.html

一些归结为这样一个事实,即einsum是新的,并且可能试图更好地缓存对齐和其他内存访问问题,而许多旧的numpy函数专注于通过大量优化的易于移植的实现。 不过,我只是在猜测。


然而,你所做的一些并不完全是“苹果对苹果”的比较。

除了@Jamie已经说过的之外, sum为数组使用了一个更合适的累加器

例如, sum对于检查输入的类型和使用适当的累加器更加小心。 例如,请考虑以下几点:

In [1]: x = 255 * np.ones(100, dtype=np.uint8)

In [2]: x
Out[2]:
array([255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255], dtype=uint8)

请注意, sum是正确的:

In [3]: x.sum()
Out[3]: 25500

虽然einsum会给出错误的结果:

In [4]: np.einsum('i->', x)
Out[4]: 156

但是如果我们使用更少的dtype ,我们仍然会得到您所期望的结果:

In [5]: y = 255 * np.ones(100)

In [6]: np.einsum('i->', y)
Out[6]: 25500.0

现在numpy 1.8发布了,根据文档,所有的ufuncs都应该使用SSE2,我想仔细检查Seberg对SSE2的评论是否有效。

为了执行测试,创建了一个新的python 2.7安装 - numpy 1.7和1.8在运行Ubuntu的AMD opteron内核上使用标准选项使用icc编译。

这是1.8升级之前和之后的测试运行:

import numpy as np
import timeit

arr_1D=np.arange(5000,dtype=np.double)
arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500)
arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500)

print 'Summation test:'
print timeit.timeit('np.sum(arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print timeit.timeit('np.einsum("ijk->", arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print '----------------------n'


print 'Power test:'
print timeit.timeit('arr_3D*arr_3D*arr_3D',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print timeit.timeit('np.einsum("ijk,ijk,ijk->ijk", arr_3D, arr_3D, arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print '----------------------n'


print 'Outer test:'
print timeit.timeit('np.outer(arr_1D, arr_1D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print timeit.timeit('np.einsum("i,k->ik", arr_1D, arr_1D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print '----------------------n'


print 'Einsum test:'
print timeit.timeit('np.sum(arr_2D*arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print timeit.timeit('np.einsum("ij,oij->", arr_2D, arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print '----------------------n'

Numpy 1.7.1:

Summation test:
0.172988510132
0.0934836149216
----------------------

Power test:
1.93524689674
0.839519000053
----------------------

Outer test:
0.130380821228
0.121401786804
----------------------

Einsum test:
0.979052495956
0.126066613197

Numpy 1.8:

Summation test:
0.116551589966
0.0920487880707
----------------------

Power test:
1.23683619499
0.815982818604
----------------------

Outer test:
0.131808176041
0.127472200394
----------------------

Einsum test:
0.781750011444
0.129271841049

我认为这是相当确定的,SSE在时间差异中起着重要作用,应该指出的是,重复这些测试的时间仅为0.003秒。 剩余的差异应该包括在这个问题的其他答案中。


我认为这些时间表明了发生了什么事情:

a = np.arange(1000, dtype=np.double)
%timeit np.einsum('i->', a)
100000 loops, best of 3: 3.32 us per loop
%timeit np.sum(a)
100000 loops, best of 3: 6.84 us per loop

a = np.arange(10000, dtype=np.double)
%timeit np.einsum('i->', a)
100000 loops, best of 3: 12.6 us per loop
%timeit np.sum(a)
100000 loops, best of 3: 16.5 us per loop

a = np.arange(100000, dtype=np.double)
%timeit np.einsum('i->', a)
10000 loops, best of 3: 103 us per loop
%timeit np.sum(a)
10000 loops, best of 3: 109 us per loop

所以,你打电话的时候基本上有一个几乎不变的3US开销np.sum超过np.einsum ,所以他们基本上跑得快,而是一个需要一点时间来走了。 为什么会这样? 我的资金如下:

a = np.arange(1000, dtype=object)
%timeit np.einsum('i->', a)
Traceback (most recent call last):
...
TypeError: invalid data type for einsum
%timeit np.sum(a)
10000 loops, best of 3: 20.3 us per loop

不知道到底发生了什么,但似乎np.einsum正在跳过一些检查来提取特定于类型的函数来执行乘法和加法运算,并且仅针对标准C类型直接使用*+


多维情况并没有不同:

n = 10; a = np.arange(n**3, dtype=np.double).reshape(n, n, n)
%timeit np.einsum('ijk->', a)
100000 loops, best of 3: 3.79 us per loop
%timeit np.sum(a)
100000 loops, best of 3: 7.33 us per loop

n = 100; a = np.arange(n**3, dtype=np.double).reshape(n, n, n)
%timeit np.einsum('ijk->', a)
1000 loops, best of 3: 1.2 ms per loop
%timeit np.sum(a)
1000 loops, best of 3: 1.23 ms per loop

所以一个大部分是不变的开销,一旦他们接近它,不是一个更快的运行。

链接地址: http://www.djcxy.com/p/85695.html

上一篇: Why is numpy's einsum faster than numpy's built in functions?

下一篇: Speed of calculating powers (in python)