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

最近在寫程式碼, 編一個 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-type
row_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 = 9999

def 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.000040
print((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.000063
np.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.000003
a[:, 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] = 0
print(a)
"""
[[1 2]
 [3 4]
 [5 6]]
"""

a_view[1,1] = 0
print(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 = 9999

def f1(a):                 
    for _ in range(N):
        a *= 2           # same as a[:] *= 2

def f2(b):
    for _ in range(N):
        b = 2*b

print('%f' % ((t1-t0)/N))     # f1: 0.000837
print('%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.001059
print('%f' % ((t2-t1)/N))    # 0.000000

 

選擇資料

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

a_view1 = a[1:2, 3:6]    # 切片 slice
a_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]]  # 用 mask
a_copy3 = a[[1,2], :]        # 雖然 1,2 的確連在一起了, 但是他們確實是 copy
a_copy4 = a[a[1,:] != 0, :]  # fancy indexing
a_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 = 99
indices = 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.000393
print('%f' % ((t2-t1)/N))    # 0.000569

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

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

mask = a[:, 0] < 0.5
def 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.028109
print('%f' % ((t2-t1)/N))    # 0.031013

 

非常有用的 out 引數

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

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

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

a += 1                 # 0.011219
np.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)

"""   
# a
array([(0,  0.), (0,  0.), (0,  0.)],
      dtype=[('foo', '<i4'), ('bar', '<f2')])

# b
   foo  bar
0    0  0.0
1    0  0.0
2    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.000003
print('%f' % ((t2-t1)/N))    # 0.000508

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

 

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

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