爲什麼numpy的einsum比numpy的內置函數更快？

Question

讓我們從三個dtype=np.double的數組開始。計時是在英特爾CPU上使用numpy 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 

Einsum似乎是至少快兩倍，爲np.inner，np.outer，np.kron，並且np.sum無論axes選擇。主要的例外是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（見change log）。我相信這是正確的答案，但不是能夠確認它。通過改變輸入數組的dtype並觀察速度差異以及不是每個人都遵守相同的定時趨勢這一事實，可以找到一些有限的證明。

Answer 1

現在numpy 1.8發佈了，根據文檔所有的ufuncs都應該使用SSE2，我想仔細檢查Seberg對SSE2的評論是否有效。

爲了執行測試，創建了一個新的python 2.7安裝 - 在運行Ubuntu的AMD opteron內核上使用標準選項，使用icc編譯numpy 1.7和1.8。

這是試運行之前和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 

我認爲這是相當確定的是，上證所在時間差異中起着重要作用，應該指出的是，重複這些測試的時機只有〜0.003秒。剩餘的差異應該包括在這個問題的其他答案中。

Answer 2

我覺得這些時間解釋這是怎麼回事：

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 

所以你基本上調用np.sum超過np.einsum當幾乎恆定3US開銷，所以他們基本上跑得快，而是一個需要一點時間來獲得去。爲什麼會這樣？我的錢是在以下方面：

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 

所以一大部分是恆定的開銷，而不是更快的運行，一旦他們着手進行。

Answer 3

首先，關於這個在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 

但是，如果我們使用限制較少我們仍然會得到您所期望的結果：

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

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