[轉]numpy性能優化

轉自：http://blog.csdn.net/pipisorry/article/details/39087583

http://blog.csdn.net/pipisorry/article/details/39087583

Introduction

NumPy提供了一個特殊的數據類型ndarray，其在向量計算上作了優化。這個對象是科學數值計算中大多數算法的核心。

相比於原生的Python，利用NumPy數組能夠得到顯著的性能加速，尤爲是當你的計算遵循單指令多數據流（SIMD）範式時。

然而，利用NumPy也有可能有意無心地寫出未優化的代碼。下面這些技巧能夠幫助你編寫高效的NumPy代碼。

避免沒必要要的數據拷貝

查看數組的內存地址

1. 查看靜默數組拷貝的第一步是在內存中找到數組的地址。下邊的函數就是作這個的：

def  id (x):      # This function returns the memory block address of an array.      return x.__array_interface__[ 'data' ][ 0 ]

2. 有時你可能須要複製一個數組，例如你須要在操做一個數組時，內存中仍然保留其原始副本。

a  = np.zeros( 10 ); aid = id (a); aid 71211328 b  = a.copy(); id (b) = = aid False

Note:python列表的淺拷貝和深拷貝-列表拷貝一節

具備相同數據地址（好比id函數的返回值）的兩個數組，共享底層數據緩衝區。然而，共享底層數據緩衝區的數組，只有當它們具備相同的偏移量（意味着它們的第一個元素相同）時，才具備相同的數據地址。共享數據緩衝區，但偏移量不一樣的兩個數組，在內存地址上有細微的差異：

id (a), id (a[ 1 :]) ( 71211328 , 71211336 )

在這篇文章中，咱們將確保函數用到的數組具備相同的偏移量。

下邊是一個判斷兩個數組是否共享相同數據的更可靠的方案：

def get_data_base(arr):      "" "For a given Numpy array, finds the base array that " owns " the actual data." ""      base = arr      while is instance(base.base, np.ndarray):          base = base.base      return  base   def arrays_share_data(x, y):      return  get_data_base(x) is  get_data_base(y)   print(arrays_share_data(a,a.copy()), arrays_share_data(a,a[ 1 :])) False True

感謝Michael Droettboom指出這種更精確的方法，提出這個替代方案。

Note:a和a[1:]id雖然不一樣，可是他們是共享內存的。

就地操做和隱式拷貝操做

3. 數組計算包括就地操做（下面第一個例子：數組修改）或隱式拷貝操做（第二個例子：建立一個新的數組）。

1 2 3 4 5

a  * = 2 ; id (a) = = aid True   c  = a * 2 ; id (c) = = aid False

必定要選擇真正須要的操做類型。隱式拷貝操做很明顯很慢，以下所示：

1 2 3 4 5 6 7

% % timeit a = np.zeros( 10000000 ) a  * = 2 10  loops, best of  3 : 19.2 ms per loop   % % timeit a = np.zeros( 10000000 ) b  = a * 2 10  loops, best of  3 : 42.6 ms per loop

4. 重塑數組可能涉及到拷貝操做，也可能涉及不到。

例如，重塑一個二維矩陣不涉及拷貝操做，除非它被轉置（或更通常的非連續操做）：

1 2	a  = np.zeros(( 10 , 10 )); aid  =  id (a); aid 53423728

重塑一個數組，同時保留其順序，並不觸發拷貝操做。

1 2	b  = a.reshape(( 1 , - 1 )); id (b) = = aid True

轉置一個數組會改變其順序，因此這種重塑會觸發拷貝操做。

1 2	c  = a.T.reshape(( 1 , - 1 )); id (c) = = aid False

所以，後邊的指令比前邊的指令明顯要慢。

5. 數組的flatten和revel方法將數組變爲一個一維向量（鋪平數組）。flatten方法老是返回一個拷貝後的副本，而revel方法只有當有必要時才返回一個拷貝後的副本（因此該方法要快得多，尤爲是在大數組上進行操做時）。

1 2 3 4 5 6 7 8 9 10 11

d  = a.flatten(); id (d) = = aid False   e  = a.ravel(); id (e) = = aid True   % timeit a.flatten() 1000000  loops, best of  3 : 881 ns per loop   % timeit a.ravel() 1000000 loops, best of 3 : 294 ns per loop

廣播規則

廣播規則容許你在形狀不一樣但卻兼容的數組上進行計算。換句話說，你並不老是須要重塑或鋪平數組，使它們的形狀匹配。

廣播規則描述了具備不一樣維度和/或形狀的數組仍能夠用於計算。通常的規則是：當兩個維度相等，或其中一個爲1時，它們是兼容的。NumPy使用這個規則，從後邊的維數開始，向前推導，來比較兩個元素級數組的形狀。最小的維度在內部被自動延伸，從而匹配其餘維度，但此操做並不涉及任何內存複製。

下面的例子說明了兩個向量之間進行矢量積的兩個方法：第一個方法涉及到數組的變形操做，第二個方法涉及到廣播規則。顯然第二個方法是要快得多。

 

n =1000

a =np.arange(n)
ac =a[:, np.newaxis]
ar =a[np.newaxis, :]

%timeit np.tile(ac, (1, n))* np.tile(ar, (n,1))
100 loops, best of 3:10 ms per loop

%timeit ar* ac
100 loops, best of 3:2.36 ms per loop

[numpy教程- 通用函數ufunc- 廣播規則]

 

NumPy數組進行高效的選擇

NumPy提供了多種數組分片的方式。

數組視圖涉及到一個數組的原始數據緩衝區，但具備不一樣的偏移量，形狀和步長。NumPy只容許等步長選擇（即線性分隔索引）。

NumPy還提供沿一個軸進行任意選擇的特定功能。

最後，花式索引（fancy indexing）是最通常的選擇方法，但正如咱們將要在文章中看到的那樣，它同時也是最慢的。

1. 建立一個具備不少行的數組。咱們將沿第一維選擇該數組的分片。

n, d  = 100000 , 100 a  = np.random.random_sample((n, d)); aid = id (a)

數組視圖和花式索引

2. 每10行選擇一行，這裏用到了兩個不一樣的方法（數組視圖和花式索引）。

b1  = a[:: 10 ] b2  = a[np.arange( 0 , n, 10 )] np.array_equal(b1, b2) True

3. 數組視圖指向原始數據緩衝區，而花式索引產生一個拷貝副本。

id (b1) = = aid, id (b2) = = aid ( True , False )

兩個方法的執行效率，花式索引慢好幾個數量級，由於它要複製一個大數組。

替代花式索引：索引列表

當須要沿一個維度進行非等步長選擇時，數組視圖就無能爲力了。

然而，替代花式索引的方法在這種狀況下依然存在。給定一個索引列表，NumPy的函數能夠沿一個軸執行選擇操做。

i  = np.arange( 0 , n, 10 )   b1  = a[i] b2  = np.take(a, i, axis = 0 )   np.array_equal(b1, b2) True

第二個方法更快一點：

1 2 3 4 5

% timeit a[i] 100  loops, best of  3 : 13 ms per loop   % timeit np.take(a, i, axis = 0 ) 100  loops, best of  3 : 4.87 ms per loop

替代花式索引：布爾掩碼

當沿一個軸進行選擇的索引是經過一個布爾掩碼向量指定時，compress函數能夠做爲花式索引的替代方案。

1	i  = np.random.random_sample(n) < . 5

可使用花式索引或者np.compress函數進行選擇。

b1  = a[i] b2  = np.compress(i, a, axis = 0 )   np.array_equal(b1, b2) True   % timeit a[i] 10  loops, best of  3 : 59.8 ms per loop   % timeit np.compress(i, a, axis = 0 ) 10  loops, best of  3 : 24.1 ms per loop

第二個方法一樣比花式索引快得多。

花式索引是進行數組任意選擇的最通常方法。然而，每每會存在更有效、更快的方法，應儘量首選那些方法。

當進行等步長選擇時應該使用數組視圖，但須要注意這樣一個事實：視圖涉及到原始數據緩衝區。

 

 

爲何NumPy數組如此高效？

一個NumPy數組基本上是由元數據（維數、形狀、數據類型等）和實際數據構成。數據存儲在一個均勻連續的內存塊中，該內存在系統內存（隨機存取存儲器，或RAM）的一個特定地址處，被稱爲數據緩衝區。這是和list等純Python結構的主要區別，list的元素在系統內存中是分散存儲的。這是使NumPy數組如此高效的決定性因素。

爲何這會如此重要？主要緣由是：

1. 低級語言好比C，能夠很高效的實現數組計算（NumPy的很大一部分其實是用C編寫）。例如，知道了內存塊地址和數據類型，數組計算只是簡單遍歷其中全部的元素。但在Python中使用list實現，會有很大的開銷。

2. 內存訪問模式中的空間位置訪問會產生顯著地性能提升，尤爲要感謝CPU緩存。事實上，緩存將字節塊從RAM加載到CPU寄存器。而後相鄰元素就能高效地被加載了（順序位置，或引用位置）。

3. 數據元素連續地存儲在內存中，因此NumPy能夠利用現代CPU的矢量化指令，像英特爾的SSE和AVX，AMD的XOP等。例如，爲了做爲CPU指令實現的矢量化算術計算，能夠加載在128,256或512位寄存器中的多個連續的浮點數。

4. NumPy能夠經過Intel Math Kernel Library (MKL)與高度優化的線性代數庫相連，好比BLAS和LAPACK。NumPy中一些特定的矩陣計算也多是多線程，充分利用了現代多核處理器的優點。

總之，將數據存儲在一個連續的內存塊中，根據內存訪問模式，CPU緩存和矢量化指令，能夠確保以最佳方式使用現代CPU的體系結構。

就地操做和隱式拷貝操做之間的區別

讓咱們解釋一下技巧3。相似於a *= 2這樣的表達式對應一個就地操做，即數組的全部元素值被乘以2。相比之下，a = a*2意味着建立了一個包含a*2結果值的新數組，變量a此時指向這個新數組。舊數組變爲了無引用的，將被垃圾回收器刪除。第一種狀況中沒有發生內存分配，相反，第二種狀況中發生了內存分配。

更通常的狀況，相似於a[i:j]這樣的表達式是數組某些部分的視圖：它們指向包含數據的內存緩衝區。利用就地操做改變它們，會改變原始數據。所以，a[:] = a * 2的結果是一個就地操做，和a = a * 2不同。

知道NumPy的這種細節能夠幫助你解決一些錯誤（例如數組由於在一個視圖上的一個操做，被無心中修改），並能經過減小沒必要要的副本數量，優化代碼的速度和內存消耗。

爲何有些數組不進行拷貝操做，就不能被重塑？

一個轉置的二維矩陣不依靠拷貝就沒法進行鋪平。一個二維矩陣包含的元素經過兩個數字（行和列）進行索引，但它在內部是做爲一個一維連續內存塊存儲的，可以使用一個數字訪問。

有多個在一維內存塊中存儲矩陣元素的方法：咱們能夠先放第一行的元素，而後第二行，以此類推，或者先放第一列的元素，而後第二列，以此類推。第一種方法叫作行優先排序，然後一種方法稱爲列優先排序。這兩種方法之間的選擇只是一個內部約定問題：NumPy使用行優先排序，相似於C，而不一樣於FORTRAN。

更通常的狀況，NumPy使用步長的概念進行多維索引和元素的底層序列（一維）內存位置之間的轉換。array[i1, i2]和內部數據的相關字節地址之間的具體映射關係爲：

1	offset = array.strides[0] * i1 + array.strides[1] * i2

重塑一個數組時，NumPy會盡量經過修改步長屬性來避免拷貝。例如，當轉置一個矩陣時，步長的順序被翻轉，但底層數據仍然是相同的。然而，僅簡單地依靠修改步長沒法完成鋪平一個轉置數組的操做（嘗試下！），因此須要一個副本。

Recipe 4.6（NumPy中使用步長技巧）包含步長方面更普遍的討論。同時，Recipe4.7（使用步長技巧實現一個高效的移動平均算法）展現瞭如何使用步伐加快特定數組計算。

from:http://blog.csdn.net/pipisorry/article/details/39087583

ref:Getting the Best Performance out of NumPy

IPython Cookbook

如何寫出比 MATLAB 更快的矩陣運算程序？

Numpy使用MKL庫提高計算性能