為什麼用 Numpy 還是慢, 你用對了嗎?

01-26

最近在寫代碼, 編一個 Python 模擬器, 做 simulation, 好不容易用傳說中 Python 里速度最快的計算模塊 "Numpy" 的寫好了, 結果運行起來, 出奇的慢! 因為一次simulation要一個小時, 要不停測試, 所以自己受不了了.. 首先, 我的腦海中的問題, 漸漸浮現出來.

我知道 Pandas 要比 Numpy 慢, 所以我盡量避免用 Pandas. 但是 Numpy (速度怪獸), 為什麼還是這麼慢?

帶有寫代碼潔癖的我好好給 google 了一番. 第一個出現在我眼前的就是這個文章, Getting the Best Performance out of NumPy. 所以我也將自己從這個文章中學到的訣竅分享給大家, 並補充一些內容.

為什麼用 Numpy?

我們都知道, Python 是慢的, 簡單來說, 因為 Python 執行你代碼的時候會執行很多複雜的 "check" 功能, 比如當你賦值

b=1; a=b/0.5

這個運算看似簡單, 但是在計算機內部, b 首先要從一個整數 integer 轉換成浮點數 float, 才能進行後面的 `b/0.5`, 因為得到的要是一個小數. 還有很多其他的原因和詳細說明 (比如 Python 對內存的調用) 在這裡能夠找到: Why Python is Slow: Looking Under the Hood

提到 Numpy, 它就是一個 Python 的救星. 能把簡單好用的 Python 和高性能的 C 語言合併在一起. 當你調用 Numpy 功能的時候, 他其實調用了很多 C 語言而不是純 Python. 這就是為什麼大家都愛用 Numpy 的原因.

創建 Numpy Array 的結構

其實 Numpy 就是 C 的邏輯, 創建存儲容器 "Array" 的時候是尋找內存上的一連串區域來存放, 而 Python 存放的時候則是不連續的區域, 這使得 Python 在索引這個容器里的數據時不是那麼有效率. Numpy 只需要再這塊固定的連續區域前後走走就能不費吹灰之力拿到數據. 下圖是來自 Why Python is Slow: Looking Under the Hood, 他很好的解釋了這一切.

在運用 Numpy 的時候, 我們通常不是用一個一維 Array 來存放數據, 而是用二維或者三維的塊來存放 (說出了學機器學習的朋友們的心聲~).

因為 Numpy 快速的矩陣相乘運算, 能將乘法運算分配到計算機中的多個核, 讓運算並行. 這年頭, 我們什麼都想多線程/多進程 (再次說出了機器學習同學們的心聲~). 這也是 Numpy 為什麼受人喜歡的一個原因. 這種並行運算大大加速了運算速度.

那麼對於這種天天要用到的2D/3D Array, 我們通常都不會想著他是怎麼來的. 因為按照我們正常人的想法, 這矩陣就是矩陣, 沒什麼深度的東西呀. 不過這可不然! 要不然我也不會寫這篇分享了. 重點來了, 不管是1D/2D/3D 的 Array, 從根本上, 它都是一個 1D array!

這篇 Blog的圖顯示. 在我們看來的 2D Array, 如果追溯到計算機內存里, 它其實是儲存在一個連續空間上的. 而對於這個連續空間, 我們如果創建 Array 的方式不同, 在這個連續空間上的排列順序也有不同. 這將影響之後所有的事情! 我們後面會用 Python 進行運算時間測試.

在 Numpy 中, 創建 2D Array 的默認方式是 "C-type" 以 row 為主在內存中排列, 而如果是 "Fortran" 的方式創建的, 就是以 column 為主在內存中排列.

col_major = np.zeros((10,10), order="C") # C-typerow_major = np.zeros((10,10), order="F") # Fortran

在 axis 上的動作

當你的計算中涉及合併矩陣, 不同形式的矩陣創建方式會給你不同的時間效果. 因為在 Numpy 中的矩陣合併等, 都是發生在一維空間里, ! 不是我們想像的二維空間中!

a = np.zeros((200, 200), order="C")b = np.zeros((200, 200), order="F")N = 9999def f1(a): for _ in range(N): np.concatenate((a, a), axis=0)def f2(b): for _ in range(N): np.concatenate((b, b), axis=0)t0 = time.time()f1(a)t1 = time.time()f2(b)t2 = time.time()print((t1-t0)/N) # 0.000040print((t2-t1)/N) # 0.000070

從上面的那張圖, 可以想到, row 為主的存儲方式, 如果在 row 的方向上合併矩陣, 將會更快. 因為只要我們將思維放在 1D array 那, 直接再加一個 row 放在1D array 後面就好了, 所以在上面的測試中, f1 速度要更快. 但是在以 column 為主的系統中, 往 1D array 後面加 row 的規則變複雜了, 消耗的時間也變長. 如果以 axis=1 的方式合併, "F" 方式的 f2 將會比 "C" 方式的 f1 更好.

還有一個要提的事情, 為了圖方便, 有時候我會直接使用 `np.stack` 來代替 `np.concatenate`, 因為這樣可以少寫一點代碼, 不過使用上面的形式, 通過上面的測試發現是這樣. 所以之後為了速度, 我推薦還是盡量使用 `np.concatenate`.

np.vstack((a,a)) # 0.000063np.concatenate((a,a), axis=0) # 0.000040

或者有時候在某個 axis 上進行操作, 比如對上面用 "C-type" 創建的 a 矩陣選點:

indices = np.random.randint(0, 100, size=10, dtype=np.int32)a[indices, :] # 0.000003a[:, indices] # 0.000006

因為 a 是用 row 為主的形式儲存, 所以在 row 上面選數據要比在 column 上選快很多! 對於其他的 axis 的操作, 結果也類似. 所以你現在懂了吧, 看自己要在哪個 axis 上動的手腳多, 然後再創建合適於自己的矩陣形式 ("C-type"/"Fortran").

copy慢 view快

在 Numpy 中, 有兩個很重要的概念, copy 和 view. copy 顧名思義, 會將數據 copy 出來存放在內存中另一個地方, 而 view 則是不 copy 數據, 直接取源數據的索引部分. 下圖來自 Understanding SettingwithCopyWarning in pandas

上面說的是什麼意思呢? 我們直接看代碼.

a = np.arange(1, 7).reshape((3,2))a_view = a[:2]a_copy = a[:2].copy()a_copy[1,1] = 0print(a)"""[[1 2] [3 4] [5 6]]"""a_view[1,1] = 0print(a)"""[[1 2] [3 0] [5 6]]"""

簡單說, a_view 的東西全部都是 a 的東西, 動 a_view 的任何地方, a 都會被動到, 因為他們在內存中的位置是一模一樣的, 本質上就是自己. 而 a_copy 則是將 a copy 了一份, 然後把 a_copy 放在內存中的另外的地方, 這樣改變 a_copy, a 是不會被改變的.

那為什麼要提這點呢? 因為 view 不會複製東西, 速度快! 我們來測試一下速度. 下面的例子中 `a*=2` 就是將這個 view 給賦值了, 和 `a[:] *= 2` 一個意思, 從頭到尾沒有創建新的東西. 而 `b = 2*b` 中, 我們將 b 賦值給另外一個新建的 b.

a = np.zeros((1000, 1000))b = np.zeros((1000, 1000))N = 9999def f1(a): for _ in range(N): a *= 2 # same as a[:] *= 2def f2(b): for _ in range(N): b = 2*bprint("%f" % ((t1-t0)/N)) # f1: 0.000837print("%f" % ((t2-t1)/N)) # f2: 0.001346

對於 view 還有一點要提, 你是不是偶爾有時候要把一個矩陣展平, 用到 `np.flatten()` 或者 `np.ravel()`. 他倆是不同的! ravel 返回的是一個 view (謝謝評論中 @非易的提醒, 官方說如果用 ravel, 需要 copy 的時候才會被 copy , 我想這個時候可能是把 ravel 裡面 order 轉換的時候, 如 "C-type" -> "Fortran"), 而 flatten 返回的總是一個 copy. 現在你知道誰在拖你的後腿了吧! 下面的測試證明, 相比於 flatten, ravel 是神速.

def f1(a): for _ in range(N): a.flatten()def f2(b): for _ in range(N): b.ravel()print("%f" % ((t1-t0)/N)) # 0.001059print("%f" % ((t2-t1)/N)) # 0.000000

選擇數據

選擇數據的時候, 我們常會用到 view 或者 copy 的形式. 我們知道了, 如果能用到 view 的, 我們就盡量用 view, 避免 copy 數據. 那什麼時候會是 view 呢? 下面舉例的都是 view 的方式:

a_view1 = a[1:2, 3:6] # 切片 slicea_view2 = a[:100] # 同上a_view3 = a[::2] # 跳步a_view4 = a.ravel() # 上面提到了... # 我只能想到這些, 如果還有請大家在評論里提出

那哪些操作我們又會變成 copy 呢?

a_copy1 = a[[1,4,6], [2,4,6]] # 用 index 選a_copy2 = a[[True, True], [False, True]] # 用 maska_copy3 = a[[1,2], :] # 雖然 1,2 的確連在一起了, 但是他們確實是 copya_copy4 = a[a[1,:] != 0, :] # fancy indexinga_copy5 = a[np.isnan(a), :] # fancy indexing... # 我只能想到這些, 如果還有請大家在評論里提出

Numpy 給了我們很多很自由的方式選擇數據, 這些雖然都很方便, 但是如果你可以盡量避免這些操作, 你的速度可以飛起來.

在上面提到的 blog 裡面, 他提到了, 如果你還是喜歡這種 fancy indexing 的形式, 我們也是可以對它加點速的. 那個 blog 中指出了兩種方法

1. 使用 `np.take()`, 替代用 index 選數據的方法.

上面提到了如果用index 來選數據, 像 `a_copy1 = a[[1,4,6], [2,4,6]]`, 用 take 在大部分情況中會比這樣的 `a_copy1` 要快.

a = np.random.rand(1000000, 10)N = 99indices = np.random.randint(0, 1000000, size=10000)def f1(a): for _ in range(N): _ = np.take(a, indices, axis=0)def f2(b): for _ in range(N): _ = b[indices]print("%f" % ((t1-t0)/N)) # 0.000393print("%f" % ((t2-t1)/N)) # 0.000569

2. 使用 `np.compress()`, 替代用 mask 選數據的方法.

上面的 `a_copy2 = a[[True, True], [False, True]]` 這種就是用 TRUE, FALSE 來選擇數據的. 測試如下:

mask = a[:, 0] < 0.5def f1(a): for _ in range(N): _ = np.compress(mask, a, axis=0)def f2(b): for _ in range(N): _ = b[mask]print("%f" % ((t1-t0)/N)) # 0.028109print("%f" % ((t2-t1)/N)) # 0.031013

非常有用的 out 參數

不深入了解 numpy 的朋友, 應該會直接忽略很多功能中的這個 out 參數 (之前我從來沒用過). 不過當我深入了解了以後, 發現他非常有用! 比如下面兩個其實在功能上是沒差的, 不過運算時間上有差, 我覺得可能是 `a=a+1` 要先轉換成 `np.add()` 這種形式再運算, 所以前者要用更久一點的時間.

a = a + 1 # 0.035230a = np.add(a, 1) # 0.032738

如果是上面那樣, 我們就會觸發之前提到的 copy 原則, 這兩個被賦值的 a, 都是原來 a 的一個 copy, 並不是 a 的 view. 但是在功能裡面有一個 out 參數, 讓我們不必要重新創建一個 a. 所以下面兩個是一樣的功能, 都不會創建另一個 copy. 不過可能是上面提到的那個原因, 這裡的運算時間也有差.

a += 1 # 0.011219np.add(a, 1, out=a) # 0.008843

帶有 out 的 numpy 功能都在這裡: Universal functions. 所以只要是已經存在了一個 placeholder (比如 a), 我們就沒有必要去再創建一個, 用 out 方便又有效.

給數據一個名字

我喜歡用 pandas, 因為 pandas 能讓你給數據命名, 用名字來做 index. 在數據類型很多的時候, 名字總是比 index 好記太多了, 也好用太多了. 但是 pandas 的確比 numpy 慢. 好在我們還是有途徑可以實現用名字來索引. 這就是 structured array. 下面 a/b 的結構是一樣的, 只是一個是 numpy 一個是 pandas.

a = np.zeros(3, dtype=[("foo", np.int32), ("bar", np.float16)])b = pd.DataFrame(np.zeros((3, 2), dtype=np.int32), columns=["foo", "bar"])b["bar"] = b["bar"].astype(np.float16)""" # aarray([(0, 0.), (0, 0.), (0, 0.)], dtype=[("foo", "<i4"), ("bar", "<f2")])# b foo bar0 0 0.01 0 0.02 0 0.0"""def f1(a): for _ in range(N): a["bar"] *= a["foo"]def f2(b): for _ in range(N): b["bar"] *= b["foo"]print("%f" % ((t1-t0)/N)) # 0.000003print("%f" % ((t2-t1)/N)) # 0.000508

可以看出來, numpy 明顯比 pandas 快很多. 如果需要使用到不同數據形式, numpy 也是可以勝任的, 並且在還保持了快速的計算速度. 至於 pandas 為什麼比 numpy 慢, 因為 pandas data 裡面還有很多七七八八的數據, 記錄著這個 data 的種種其他的特徵. 這裡還有更全面的對比: Numpy Vs Pandas Performance Comparison

如果大家還有其他的小技巧或者是速度大比拼, 歡迎在下面討論. (一切為了速度~)

最後插個小廣告, 如果你對機器學習感興趣, 這裡有很多厲害的短片形式機器學習方法介紹和很多機器學習的 Python 實踐教程, 讓你可以用業餘時間秒懂機器學習: 莫煩Python 教程