NumPy是Python中眾多科學(xué)軟件包的基礎(chǔ)。它提供了一個特殊的數(shù)據(jù)類型ndarray，其在向量計(jì)算上做了優(yōu)化。這個對象是科學(xué)數(shù)值計(jì)算中大多數(shù)算法的核心。

相比于原生的Python，利用NumPy數(shù)組可以獲得顯著的性能加速，尤其是當(dāng)你的計(jì)算遵循單指令多數(shù)據(jù)流（SIMD）范式時(shí)。然而，利用NumPy也有可能有意無意地寫出未優(yōu)化的代碼。

在這篇文章中，我們將看到一些技巧，這些技巧可以幫助你編寫高效的NumPy代碼。我們首先看一下如何避免不必要的數(shù)組拷貝，以節(jié)省時(shí)間和內(nèi)存。因此，我們將需要深入NumPy的內(nèi)部。
學(xué)習(xí)避免不必要的數(shù)據(jù)拷貝

NumPy數(shù)組計(jì)算可能涉及到內(nèi)存塊之間的內(nèi)部拷貝。有時(shí)會有不必要的拷貝，此時(shí)應(yīng)該避免。相應(yīng)地這里有一些技巧，可以幫助你優(yōu)化你的代碼。

import numpy as np

查看數(shù)組的內(nèi)存地址

1. 查看靜默數(shù)組拷貝的第一步是在內(nèi)存中找到數(shù)組的地址。下邊的函數(shù)就是做這個的：
 

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

2. 有時(shí)你可能需要復(fù)制一個數(shù)組，例如你需要在操作一個數(shù)組時(shí)，內(nèi)存中仍然保留其原始副本。
 

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

具有相同數(shù)據(jù)地址（比如id函數(shù)的返回值）的兩個數(shù)組，共享底層數(shù)據(jù)緩沖區(qū)。然而，共享底層數(shù)據(jù)緩沖區(qū)的數(shù)組，只有當(dāng)它們具有相同的偏移量（意味著它們的第一個元素相同）時(shí)，才具有相同的數(shù)據(jù)地址。共享數(shù)據(jù)緩沖區(qū)，但偏移量不同的兩個數(shù)組，在內(nèi)存地址上有細(xì)微的差別，正如下邊的例子所展示的那樣：
 

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

在這篇文章中，我們將確保函數(shù)用到的數(shù)組具有相同的偏移量。下邊是一個判斷兩個數(shù)組是否共享相同數(shù)據(jù)的更可靠的方案：

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. 數(shù)組計(jì)算包括就地操作（下面第一個例子：數(shù)組修改）或隱式拷貝操作（第二個例子：創(chuàng)建一個新的數(shù)組）。
 

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

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

%%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. 重塑一個數(shù)組可能涉及到拷貝操作，也可能涉及不到。原因?qū)⒃谙旅娼忉?。例如，重塑一個二維矩陣不涉及拷貝操作，除非它被轉(zhuǎn)置（或更一般的非連續(xù)操作）：

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

重塑一個數(shù)組，同時(shí)保留其順序，并不觸發(fā)拷貝操作。
 

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

轉(zhuǎn)置一個數(shù)組會改變其順序，所以這種重塑會觸發(fā)拷貝操作。
 

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

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

5. 數(shù)組的flatten和revel方法將數(shù)組變?yōu)橐粋€一維向量（鋪平數(shù)組）。flatten方法總是返回一個拷貝后的副本，而revel方法只有當(dāng)有必要時(shí)才返回一個拷貝后的副本（所以該方法要快得多，尤其是在大數(shù)組上進(jìn)行操作時(shí)）。
 

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

廣播規(guī)則

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

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數(shù)組上進(jìn)行高效的選擇

NumPy提供了多種數(shù)組分片的方式。數(shù)組視圖涉及到一個數(shù)組的原始數(shù)據(jù)緩沖區(qū)，但具有不同的偏移量，形狀和步長。NumPy只允許等步長選擇（即線性分隔索引）。NumPy還提供沿一個軸進(jìn)行任意選擇的特定功能。最后，花式索引（fancy indexing）是最一般的選擇方法，但正如我們將要在文章中看到的那樣，它同時(shí)也是最慢的。如果可能，我們應(yīng)該選擇更快的替代方法。

1. 創(chuàng)建一個具有很多行的數(shù)組。我們將沿第一維選擇該數(shù)組的分片。
 

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

數(shù)組視圖和花式索引

2. 每10行選擇一行，這里用到了兩個不同的方法（數(shù)組視圖和花式索引）。
 

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

3. 數(shù)組視圖指向原始數(shù)據(jù)緩沖區(qū)，而花式索引產(chǎn)生一個拷貝副本。
 

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

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

%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

花式索引慢好幾個數(shù)量級，因?yàn)樗獜?fù)制一個大數(shù)組。
替代花式索引：索引列表

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

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

第二個方法更快一點(diǎn)：
 

%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. 當(dāng)沿一個軸進(jìn)行選擇的索引是通過一個布爾掩碼向量指定時(shí)，compress函數(shù)可以作為花式索引的替代方案。
 

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

可以使用花式索引或者np.compress函數(shù)進(jìn)行選擇。
 

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

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

花式索引是進(jìn)行數(shù)組任意選擇的最一般方法。然而，往往會存在更有效、更快的方法，應(yīng)盡可能首選那些方法。

當(dāng)進(jìn)行等步長選擇時(shí)應(yīng)該使用數(shù)組視圖，但需要注意這樣一個事實(shí)：視圖涉及到原始數(shù)據(jù)緩沖區(qū)。
它是如何工作的？

在本節(jié)中，我們將看到使用NumPy時(shí)底層會發(fā)生什么，從而讓我們理解該文章中的優(yōu)化技巧。
為什么NumPy數(shù)組如此高效？

一個NumPy數(shù)組基本上是由元數(shù)據(jù)（維數(shù)、形狀、數(shù)據(jù)類型等）和實(shí)際數(shù)據(jù)構(gòu)成。數(shù)據(jù)存儲在一個均勻連續(xù)的內(nèi)存塊中，該內(nèi)存在系統(tǒng)內(nèi)存（隨機(jī)存取存儲器，或RAM）的一個特定地址處，被稱為數(shù)據(jù)緩沖區(qū)。這是和list等純Python結(jié)構(gòu)的主要區(qū)別，list的元素在系統(tǒng)內(nèi)存中是分散存儲的。這是使NumPy數(shù)組如此高效的決定性因素。

為什么這會如此重要？主要原因是：

1. 低級語言比如C，可以很高效的實(shí)現(xiàn)數(shù)組計(jì)算（NumPy的很大一部分實(shí)際上是用C編寫）。例如，知道了內(nèi)存塊地址和數(shù)據(jù)類型，數(shù)組計(jì)算只是簡單遍歷其中所有的元素。但在Python中使用list實(shí)現(xiàn)，會有很大的開銷。

2. 內(nèi)存訪問模式中的空間位置訪問會產(chǎn)生顯著地性能提高，尤其要感謝CPU緩存。事實(shí)上，緩存將字節(jié)塊從RAM加載到CPU寄存器。然后相鄰元素就能高效地被加載了（順序位置，或引用位置）。

3. 數(shù)據(jù)元素連續(xù)地存儲在內(nèi)存中，所以NumPy可以利用現(xiàn)代CPU的矢量化指令，像英特爾的SSE和AVX，AMD的XOP等。例如，為了作為CPU指令實(shí)現(xiàn)的矢量化算術(shù)計(jì)算，可以加載在128,256或512位寄存器中的多個連續(xù)的浮點(diǎn)數(shù)。

此外，說一下這樣一個事實(shí)：NumPy可以通過Intel Math Kernel Library (MKL)與高度優(yōu)化的線性代數(shù)庫相連，比如BLAS和LAPACK。NumPy中一些特定的矩陣計(jì)算也可能是多線程，充分利用了現(xiàn)代多核處理器的優(yōu)勢。

總之，將數(shù)據(jù)存儲在一個連續(xù)的內(nèi)存塊中，根據(jù)內(nèi)存訪問模式，CPU緩存和矢量化指令，可以確保以最佳方式使用現(xiàn)代CPU的體系結(jié)構(gòu)。
就地操作和隱式拷貝操作之間的區(qū)別是什么？

讓我們解釋一下技巧3。類似于a *= 2這樣的表達(dá)式對應(yīng)一個就地操作，即數(shù)組的所有元素值被乘以2。相比之下，a = a*2意味著創(chuàng)建了一個包含a*2結(jié)果值的新數(shù)組，變量a此時(shí)指向這個新數(shù)組。舊數(shù)組變?yōu)榱藷o引用的，將被垃圾回收器刪除。第一種情況中沒有發(fā)生內(nèi)存分配，相反，第二種情況中發(fā)生了內(nèi)存分配。

更一般的情況，類似于a[i:j]這樣的表達(dá)式是數(shù)組某些部分的視圖：它們指向包含數(shù)據(jù)的內(nèi)存緩沖區(qū)。利用就地操作改變它們，會改變原始數(shù)據(jù)。因此，a[:] = a * 2的結(jié)果是一個就地操作，和a = a * 2不一樣。

知道NumPy的這種細(xì)節(jié)可以幫助你解決一些錯誤（例如數(shù)組因?yàn)樵谝粋€視圖上的一個操作，被無意中修改），并能通過減少不必要的副本數(shù)量，優(yōu)化代碼的速度和內(nèi)存消耗。
為什么有些數(shù)組不進(jìn)行拷貝操作，就不能被重塑？

我們在這里解釋下技巧4，一個轉(zhuǎn)置的二維矩陣不依靠拷貝就無法進(jìn)行鋪平。一個二維矩陣包含的元素通過兩個數(shù)字（行和列）進(jìn)行索引，但它在內(nèi)部是作為一個一維連續(xù)內(nèi)存塊存儲的，可使用一個數(shù)字訪問。有多個在一維內(nèi)存塊中存儲矩陣元素的方法：我們可以先放第一行的元素，然后第二行，以此類推，或者先放第一列的元素，然后第二列，以此類推。第一種方法叫做行優(yōu)先排序，而后一種方法稱為列優(yōu)先排序。這兩種方法之間的選擇只是一個內(nèi)部約定問題：NumPy使用行優(yōu)先排序，類似于C，而不同于FORTRAN。

更一般的情況，NumPy使用步長的概念進(jìn)行多維索引和元素的底層序列（一維）內(nèi)存位置之間的轉(zhuǎn)換。array[i1, i2]和內(nèi)部數(shù)據(jù)的相關(guān)字節(jié)地址之間的具體映射關(guān)系為：
 

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

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

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

內(nèi)部數(shù)組排列還可以解釋一些NumPy相似操作之間的意想不到的性能差異。作為一個小練習(xí)，你能解釋一下下邊這個例子嗎？
 

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

NumPy的廣播規(guī)則是什么？

廣播規(guī)則描述了具有不同維度和/或形狀的數(shù)組仍可以用于計(jì)算。一般的規(guī)則是：當(dāng)兩個維度相等，或其中一個為1時(shí)，它們是兼容的。NumPy使用這個規(guī)則，從后邊的維數(shù)開始，向前推導(dǎo)，來比較兩個元素級數(shù)組的形狀。最小的維度在內(nèi)部被自動延伸，從而匹配其他維度，但此操作并不涉及任何內(nèi)存復(fù)制。

最新評論