Python中最佳化NumPy包使用效能的教程

NumPy是Python中眾多科學軟體包的基礎。它提供了一個特殊的資料類型ndarray，其在向量計算上做了最佳化。這個對象是科學數值計算中大多數演算法的核心。

相比於原生的Python，利用NumPy數組可以獲得顯著的效能加速，尤其是當你的計算遵循單指令多資料流（SIMD）範式時。然而，利用NumPy也有可能有意無意地寫出未最佳化的代碼。

在這篇文章中，我們將看到一些技巧，這些技巧可以協助你編寫高效的NumPy代碼。我們首先看一下如何避免不必要的數組拷貝，以節省時間和記憶體。因此，我們將需要深入NumPy的內部。
學習避免不必要的資料拷貝

NumPy數組計算可能涉及到記憶體塊之間的內部拷貝。有時會有不必要的拷貝，此時應該避免。相應地這裡有一些技巧，可以協助你最佳化你的代碼。

import numpy as np

查看數組的記憶體位址

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); aid71211328b = a.copy(); id(b) == aidFalse

具有相同資料地址（比如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 isinstance(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指出這種更精確的方法，提出這個替代方案。
就地操作和隱式拷貝操作

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

a *= 2; id(a) == aidTrue c = a * 2; id(c) == aidFalse

一定要選擇真正需要的操作類型。隱式拷貝操作很明顯很慢，如下所示：

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

4. 重塑一個數組可能涉及到拷貝操作，也可能涉及不到。原因將在下面解釋。例如，重塑一個二維矩陣不涉及拷貝操作，除非它被轉置（或更一般的非連續操作）：

a = np.zeros((10, 10)); aid = id(a); aid53423728

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

b = a.reshape((1, -1)); id(b) == aidTrue

轉置一個數組會改變其順序，所以這種重塑會觸發拷貝操作。

c = a.T.reshape((1, -1)); id(c) == aidFalse

因此，後邊的指令比前邊的指令明顯要慢。

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

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

廣播規則

6. 廣播規則允許你在形狀不同但卻相容的數組上進行計算。換句話說，你並不總是需要重塑或鋪平數組，使它們的形狀匹配。下面的例子說明了兩個向量之間進行向量積的兩個方法：第一個方法涉及到數組的變形操作，第二個方法涉及到廣播規則。顯然第二個方法是要快得多。

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 * ac100 loops, best of 3: 2.36 ms per loop

在NumPy數組上進行高效的選擇

NumPy提供了多種數組分區的方式。數組視圖涉及到一個數組的未經處理資料緩衝區，但具有不同的位移量，形狀和步長。NumPy只允許等步長選擇（即線性分隔索引）。NumPy還提供沿一個軸進行任意選擇的特定功能。最後，花式索引（fancy indexing）是最一般的選擇方法，但正如我們將要在文章中看到的那樣，它同時也是最慢的。如果可能，我們應該選擇更快的替代方法。

1. 建立一個具有很多行的數組。我們將沿第一維選擇該數組的分區。

n, d = 100000, 100a = 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)

4. 比較一下兩個方法的執行效率。

%timeit a[::10]1000000 loops, best of 3: 804 ns per loop %timeit a[np.arange(0, n, 10)]100 loops, best of 3: 14.1 ms per loop

花式索引慢好幾個數量級，因為它要複製一個大數組。
替代花式索引：索引列表

5. 當需要沿一個維度進行非等步長選擇時，數組視圖就無能為力了。然而，替代花式索引的方法在這種情況下依然存在。給定一個索引列表，NumPy的函數可以沿一個軸執行選擇操作。

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

第二個方法更快一點：

%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

替代花式索引：布爾掩碼

6. 當沿一個軸進行選擇的索引是通過一個布爾掩碼向量指定時，compress函數可以作為花式索引的替代方案。

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數組如此高效？

一個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位寄存器中的多個連續的浮點數。

此外，說一下這樣一個事實：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的這種細節可以協助你解決一些錯誤（例如數組因為在一個視圖上的一個操作，被無意中修改），並能通過減少不必要的副本數量，最佳化代碼的速度和記憶體消耗。
為什麼有些數組不進行拷貝操作，就不能被重塑？

我們在這裡解釋下技巧4，一個轉置的二維矩陣不依靠拷貝就無法進行鋪平。一個二維矩陣包含的元素通過兩個數字（行和列）進行索引，但它在內部是作為一個一維連續記憶體Block Storage的，可使用一個數字訪問。有多個在一維記憶體塊中儲存矩陣元素的方法：我們可以先放第一行的元素，然後第二行，以此類推，或者先放第一列的元素，然後第二列，以此類推。第一種方法叫做行優先排序，而後一種方法稱為列優先排序。這兩種方法之間的選擇只是一個內部約定問題：NumPy使用行優先排序，類似於C，而不同於FORTRAN。

更一般的情況，NumPy使用步長的概念進行多維索引和元素的底層序列（一維）記憶體位置之間的轉換。array[i1, i2]和內部資料的相關位元組地址之間的具體映射關係為：

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

重塑一個數組時，NumPy會儘可能通過修改步長屬性來避免拷貝。例如，當轉置一個矩陣時，步長的順序被翻轉，但底層資料仍然是相同的。然而，僅簡單地依靠修改步長無法完成鋪平一個轉置數組的操作（嘗試下！），所以需要一個副本。

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

內部數組排列還可以解釋一些NumPy相似操作之間的意想不到的效能差異。作為一個小練習，你能解釋一下下邊這個例子嗎？

a = np.random.rand(5000, 5000)%timeit a[0,:].sum()%timeit a[:,0].sum() 100000 loops, best of 3: 9.57 μs per loop10000 loops, best of 3: 68.3 μs per loop

NumPy的廣播規則是什嗎？

廣播規則描述了具有不同維度和/或形狀的數組仍可以用於計算。一般的規則是：當兩個維度相等，或其中一個為1時，它們是相容的。NumPy使用這個規則，從後邊的維數開始，向前推導，來比較兩個元素級數組的形狀。最小的維度在內部被自動延伸，從而匹配其他維度，但此操作並不涉及任何記憶體複製。

