如何讓Python像Julia一樣快地運行
Julia 與 Python 的比較
我是否應丟棄 Python 和其他語言，使用 Julia 執行技術計算?在看到 http://julialang.org/ 上的基準測試后，人們一定會這么想。Python和其他高級語言在速度上遠遠有些落后。但是，我想到的第一個問題有所不同：Julia 團隊能否以最適合 Python 的方式編寫 Python 基準測試?
我對這種跨語言比較的觀點是，應該根據要執行的任務來定義基準測試，然后由語言專家編寫執行這些任務的最佳代碼。如果代碼全由一個語言團隊編寫，則存在其他語言未得到最佳使用的風險。
Julia 團隊有一件事做得對，那就是他們將他們使用的代碼發布到了 github 上。具體地講，Python 代碼可在此處找到。
第一眼看到該代碼，就可以證實我所害怕的偏見。該代碼是以 C 風格編寫的，在數組和列表上大量使用了循環。這不是使用 Python 的最佳方式。
我不會責怪 Julia 團隊，因為我很內疚自己也有同樣的偏見。但我受到了殘酷的教訓：付出任何代價都要避免數組或列表上的循環，因為它們確實會拖慢 Python中的速度，請參閱《Python 不是 C》。
考慮到對 C 風格的這種偏見，一個有趣的問題(至少對我而言)是，我們能否改進這些基準測試，更好地使用 Python 及其工具?
在我給出答案之前，我想說我絕不會試圖貶低 Julia。在進一步開發和改進后，Julia 無疑是一種值得關注的語言。我只是想分析 Python方面的事情。實際上，我正在以此為借口來探索各種可用于讓代碼更快運行的 Python 工具。
在下面的內容中，我使用 Docker 鏡像在 Jupyter Notebook 中使用 Python 3.4.3，其中已安裝了所有的 Python 科學工具組合。我還會通過Windows 機器上的 Python 2.7.10，使用 Anaconda 來運行代碼。計時是對 Python 3.4.3 執行的。包含下面的所有基準測試的完整代碼的 Notebook 可在此處(https://www.ibm.com/developerworks/community/blogs/jfp/resource/julia_python.zip)找到。
鑒于各種社交媒體上的評論，我添加了這樣一句話：我沒有在這里使用 Python 的替代性實現。我沒有編寫任何 C代碼：如果您不信，可試試尋找分號。本文中使用的所有工具都是 Anaconda 或其他發行版中提供的標準的 Cython 實現。下面的所有代碼都在單個 Notebook中運行。
我嘗試過使用來自 github 的 Julia 微性能文件，但不能使用 Julia 0.4.2 原封不動地運行它。我必須編輯它并將 @timeit 替換為@time，它才能運行。在對它們計時之前，我還必須添加對計時函數的調用，否則編譯時間也將包含在內。我使用的文件位于此處。我在用于運行 Python 的同一個機器上使用 Julia 命令行接口運行它。
計時代碼
Julia 團隊使用的第一項基準測試是 Fibonacci 函數的一段簡單編碼。
def fib(n):
if n<2:
return n
return fib(n-1)+fib(n-2)
此函數的值隨 n 的增加而快速增加，例如：
fib(100) = 354224848179261915075
可以注意到，Python 任意精度 (arbitrary precision) 很方便。在 C 等語言中編寫相同的函數需要花一些編碼工作來避免整數溢出。在 Julia中，需要使用 BigInt 類型。
所有 Julia 基準測試都與運行時間有關。這是 Julia 中使用和不使用 BigInt 的計時：
0.000080 seconds (149 allocations:10.167 KB)
0.012717 seconds (262.69 k allocations:4.342 MB)
在 Python Notebook 中獲得運行時間的一種方式是使用神奇的 %timeit。例如，在一個新單元中鍵入：
%timeit fib(20)
執行它會獲得輸出：
100 loops, best of 3:3.33 ms per loop
這意味著計時器執行了以下操作：
- 運行 fib(20) 100 次，存儲總運行時間
- 運行 fib(20) 100 次，存儲總運行時間
- 運行 fib(20) 100 次，存儲總運行時間
- 從 3 次運行中獲取最小的運行時間，將它除以 100，然后輸出結果，該結果就是 fib(20) 的最佳運行時間
這些循環的大小(100 次和 3 次)會由計時器自動調整。可能會根據被計時的代碼的運行速度來更改循環大小。
Python 計時與使用了 BigInt 時的 Julia 計時相比出色得多：3 毫秒與 12 毫秒。在使用任意精度時，Python 的速度是 Julia 的 4倍。
但是，Python 比 Julia 默認的 64 位整數要慢。我們看看如何在 Python 中強制使用 64 位整數。
使用 Cython 編譯
一種編譯方式是使用 Cython 編譯器。這個編譯器是使用 Python
編寫的。它可以通過以下命令安裝：
pip install Cython
如果使用 Anaconda，安裝會有所不同。因為安裝有點復雜，所以我編寫了一篇相關的博客文章：將 Cython For Anaconda 安裝在 Windows 上
安裝后，我們使用神奇的 %load_ext 將 Cython 加載到 Notebook 中：
%load_ext Cython
然后就可以在我們的 Notebook 中編譯代碼。我們只需要將想要編譯的代碼放在一個單元中，包括所需的導入語句，使用神奇的 %%cython 啟動該單元：
%%cython
def fib_cython(n):
if n<2:
return n
return fib_cython(n-1)+fib_cython(n-2)
執行該單元會無縫地編譯這段代碼。我們為該函數使用一個稍微不同的名稱，以反映出它是使用 Cython編譯的。當然，一般不需要這么做。我們可以將之前的函數替換為相同名稱的已編譯函數。
對它計時會得到：
1000 loops, best of 3:1.22 ms per loop
哇，幾乎比最初的 Python 代碼快 3 倍!我們現在比使用 BigInt 的 Julia 快 100 倍。
我們還可以嘗試靜態類型。使用關鍵字 cpdef 而不是 def 來聲明該函數。它使我們能夠使用相應的 C 類型來鍵入函數的參數。我們的代碼變成了：
%%cython
cpdef long fib_cython_type(long n):
if n<2:
return n
return fib_cython_type(n-1)+fib_cython_type(n-2)
執行該單元后，對它計時會得到：
10000 loops, best of 3:36 μs per loop
太棒了，我們現在只花費了 36 微秒，比最初的基準測試快約 100 倍!這與 Julia 所花的 80 毫秒相比更出色。
有人可能會說，靜態類型違背了 Python的用途。一般來講，我比較同意這種說法，我們稍后將查看一種在不犧牲性能的情況下避免這種情形的方法。但我并不認為這是一個問題。Fibonacci函數必須使用整數來調用。我們在靜態類型中失去的是 Python 所提供的任意精度。對于 Fibonacci，使用 C 類型 long
會限制輸入參數的大小，因為太大的參數會導致整數溢出。
請注意，Julia 計算也是使用 64 位整數執行的，因此將我們的靜態類型版本與 Julia 的對比是公平的。
緩存計算
我們在保留 Python 任意精度的情況下能做得更好。fib 函數重復執行同一種計算許多次。例如，fib(20) 將調用 fib(19) 和fib(18)。fib(19) 將調用 fib(18) 和 fib(17)。結果 fib(18) 被調用了兩次。簡單分析表明，fib(17) 將被調用 3 次，fib(16) 將被調用 5 次，等等。
在 Python 3 中，我們可以使用 functools 標準庫來避免這些重復的計算。
from functools import lru_cache as cache
@cache(maxsize=None)
def fib_cache(n):
if n<2:
return n
return fib_cache(n-1)+fib_cache(n-2)
對此函數計時會得到：
1000000 loops, best of 3:910 ns per loop
速度又增加了 40 倍，比最初的 Python 代碼快約 3,600 倍!考慮到我們僅向遞歸函數添加了一條注釋，此結果非常令人難忘。
Python 2.7 中沒有提供這種自動緩存。我們需要顯式地轉換代碼，才能避免這種情況下的重復計算。
def fib_seq(n):
if n < 2:
return n
a,b = 1,0
for i in range(n-1):
a,b = a+b,a
return a
請注意，此代碼使用了 Python 同時分配兩個局部變量的能力。對它計時會得到：
1000000 loops, best of 3:1.77 μs per loop
我們又快了 20 倍!讓我們在使用和不使用靜態類型的情況下編譯我們的函數。請注意，我們使用了 cdef 關鍵字來鍵入局部變量。
%%cython
def fib_seq_cython(n):
if n < 2:
return n
a,b = 1,0
for i in range(n-1):
a,b = a+b,a
return a
cpdef long fib_seq_cython_type(long n):
if n < 2:
return n
cdef long a,b
a,b = 1,0
for i in range(n-1):
a,b = a+b,b
return a
我們可在一個單元中對兩個版本計時：
%timeit fib_seq_cython(20)
%timeit fib_seq_cython_type(20)
結果為：
1000000 loops, best of 3:953 ns per loop
10000000 loops, best of 3:51.9 ns per loop
靜態類型代碼現在花費的時間為 51.9 納秒，比最初的基準測試快約 60,000(六萬)倍。
如果我們想計算任意輸入的 Fibonacci 數，我們應堅持使用無類型版本，該版本的運行速度快 3,500 倍。還不錯，對吧?
使用 Numba 編譯
讓我們使用另一個名為 Numba 的工具。它是針對部分 Python 版本的一個即時
(jit) 編譯器。它不是對所有 Python 版本都適用，但在適用的情況下，它會帶來奇跡。
安裝它可能很麻煩。推薦使用像 Anaconda 這樣的 Python 發行版或一個已安裝了 Numba 的 Docker 鏡像。完成安裝后，我們導入它的 jit 編譯器：
from numba import jit
它的使用非常簡單。我們僅需要向想要編譯的函數添加一點修飾。我們的代碼變成了：
@jit
def fib_seq_numba(n):
if n < 2:
return n
(a,b) = (1,0)
for i in range(n-1):
(a,b) = (a+b,a)
return a
對它計時會得到：
1000000 loops, best of 3:225 ns per loop
比無類型的 Cython 代碼更快，比最初的 Python 代碼快約 16,000 倍!
使用 Numpy
我們現在來看看第二項基準測試。它是快速排序算法的實現。Julia 團隊使用了以下 Python 代碼：
def qsort_kernel(a, lo, hi):
i = lo
j = hi
while i < hi:
pivot = a[(lo+hi) // 2]
while i <= j:
while a[i] < pivot:
i += 1
while a[j] > pivot:
j -= 1
if i <= j:
a[i], a[j] = a[j], a[i]
i += 1
j -= 1
if lo < j:
qsort_kernel(a, lo, j)
lo = i
j = hi
return a
我將他們的基準測試代碼包裝在一個函數中：
import random
def benchmark_qsort():
lst = [ random.random() for i in range(1,5000) ]
qsort_kernel(lst, 0, len(lst)-1)
對它計時會得到：
100 loops, best of 3:18.3 ms per loop
上述代碼與 C 代碼非常相似。Cython 應該能很好地處理它。除了使用 Cython 和靜態類型之外，讓我們使用 Numpy
數組代替列表。在數組大小較大時，比如數千個或更多元素，Numpy 數組確實比Python 列表更快。
安裝 Numpy 可能會花一些時間，推薦使用 Anaconda 或一個已安裝了 Python 科學工具組合的 Docker 鏡像。
在使用 Cython 時，需要將 Numpy 導入到應用了 Cython 的單元中。在使用 C 類型時，還必須使用 cimport 將它作為 C 模塊導入。Numpy數組使用一種表示數組元素類型和數組維數(一維、二維等)的特殊語法來聲明。
%%cython
import numpy as np
cimport numpy as np
cpdef np.ndarray[double, ndim=1] \
qsort_kernel_cython_numpy_type(np.ndarray[double, ndim=1] a, \
long lo, \
long hi):
cdef:
long i, j
double pivot
i = lo
j = hi
while i < hi:
pivot = a[(lo+hi) // 2]
while i <= j:
while a[i] < pivot:
i += 1
while a[j] > pivot:
j -= 1
if i <= j:
a[i], a[j] = a[j], a[i]
i += 1
j -= 1
if lo < j:
qsort_kernel_cython_numpy_type(a, lo, j)
lo = i
j = hi
return a
cpdef benchmark_qsort_numpy_cython():
lst = np.random.rand(5000)
qsort_kernel_cython_numpy_type(lst, 0, len(lst)-1)
對 benchmark_qsort_numpy_cython() 函數計時會得到：
1000 loops, best of 3:1.32 ms per loop
我們比最初的基準測試快了約 15 倍，但這仍然不是使用 Python 的最佳方法。最佳方法是使用 Numpy 內置的 sort()函數。它的默認行為是使用快速排序算法。對此代碼計時：
def benchmark_sort_numpy():
lst = np.random.rand(5000)
np.sort(lst)
會得到：
1000 loops, best of 3:350 μs per loop
我們現在比最初的基準測試快 52 倍!Julia 在該基準測試上花費了 419 微秒，因此編譯的 Python 快 20%。
我知道，一些讀者會說我不會進行同類比較。我不同意。請記住，我們現在的任務是使用主機語言以最佳的方式排序輸入數組。在這種情況下，最佳方法是使用一個內置的函數。
剖析代碼
我們現在來看看第三個示例，計算 Mandelbrodt 集。Julia 團隊使用了這段 Python 代碼：
def mandel(z):
maxiter = 80
c = z
for n in range(maxiter):
if abs(z) > 2:
return n
z = z*z + c
return maxiter
def mandelperf():
r1 = np.linspace(-2.0, 0.5, 26)
r2 = np.linspace(-1.0, 1.0, 21)
return [mandel(complex(r, i)) for r in r1 for i in r2]
assert sum(mandelperf()) == 14791
最后一行是一次合理性檢查。對 mandelperf() 函數計時會得到：
100 loops, best of 3:4.62 ms per loop
使用 Cython 會得到：
100 loops, best of 3:1.94 ms per loop
還不錯，但我們可以使用 Numba 做得更好。不幸的是，Numba 還不會編譯列表推導式 (list
comprehension)。因此，我們不能將它應用到第二個函數，但我們可以將它應用到第一個函數。我們的代碼類似以下代碼。
@jit
def mandel_numba(z):
maxiter = 80
c = z
for n in range(maxiter):
if abs(z) > 2:
return n
z = z*z + c
return maxiter
def mandelperf_numba():
r1 = np.linspace(-2.0, 0.5, 26)
r2 = np.linspace(-1.0, 1.0, 21)
return [mandel_numba(complex(r, i)) for r in r1 for i in r2]
對它計時會得到：
1000 loops, best of 3:503 μs per loop
還不錯，比 Cython 快 4 倍，比最初的 Python 代碼快 9 倍!
我們還能做得更好嗎?要知道是否能做得更好，一種方式是剖析代碼。內置的 %prun 剖析器在這里不夠精確，我們必須使用一個稱為 line_profiler 的更好的剖析器。它可以通過pip 進行安裝：
pip install line_profiler
安裝后，我們需要加載它：
%load_ext line_profiler
然后使用一個神奇的命令剖析該函數：
%lprun -s -f mandelperf_numba mandelperf_numba()
它在一個彈出窗口中輸出以下信息。
Timer unit:1e-06 s
Total time:0.003666 s
File:
Function: mandelperf_numba at line 11
Line # Hits Time Per Hit % Time Line Contents
==============================================================
11 def mandelperf_numba():
12 1 1994 1994.0 54.4 r1 = np.linspace(-2.0, 0.5, 26)
13 1 267 267.0 7.3 r2 = np.linspace(-1.0, 1.0, 21)
14 1 1405 1405.0 38.3 return [mandel_numba(complex(r, i)) for r in r1 for i in r2]
我們看到，大部分時間都花費在了 mandelperf_numba() 函數的第一行和最后一行上。最后一行有點復雜，讓我們將它分為兩部分來再次剖析：
def mandelperf_numba():
r1 = np.linspace(-2.0, 0.5, 26)
r2 = np.linspace(-1.0, 1.0, 21)
c3 = [complex(r, i) for r in r1 for i in r2]
return [mandel_numba(c) for c in c3]
剖析器輸出變成：
Timer unit:1e-06 s
Total time:0.002002 s
File:
Function: mandelperf_numba at line 11
Line # Hits Time Per Hit % Time Line Contents
==============================================================
11 def mandelperf_numba():
12 1 678 678.0 33.9 r1 = np.linspace(-2.0, 0.5, 26)
13 1 235 235.0 11.7 r2 = np.linspace(-1.0, 1.0, 21)
14 1 617 617.0 30.8 c3 = [complex(r, i) for r in r1 for i in r2]
15 1 472 472.0 23.6 return [mandel_numba(c) for c in c3]
我們可以看到，對函數 mandel_numba() 的調用僅花費了總時間的 1/4。剩余時間花在 mandelperf_numba()
函數上。花時間優化它是值得的。
再次使用 Numpy
使用 Cython 在這里沒有太大幫助，而且 Numba 不適用。擺脫此困境的一種方法是再次使用 Numpy。我們將以下代碼替換為生成等效結果的 Numpy
代碼。
return [mandel_numba(complex(r, i)) for r in r1 for i in r2]
此代碼構建了所謂的二維網格。它計算由 r1 和 r2 提供坐標的點的復數表示。點 Pij 的坐標為 r1[i] 和 r2[j]。Pij 通過復數 r1[i] +
1j*r2[j] 進行表示，其中特殊常量 1j 表示單個虛數 i。
我們可以直接編寫此計算的代碼：
@jit
def mandelperf_numba_mesh():
width = 26
height = 21
r1 = np.linspace(-2.0, 0.5, width)
r2 = np.linspace(-1.0, 1.0, height)
mandel_set = np.zeros((width,height), dtype=int)
for i in range(width):
for j in range(height):
mandel_set[i,j] = mandel_numba(r1[i] + 1j*r2[j])
return mandel_set
請注意，我將返回值更改為了一個二維整數數組。如果要顯示結果，該結果與我們需要的結果更接近。
對它計時會得到：
10000 loops, best of 3:140 μs per loop
我們比最初的 Python 代碼快約 33 倍!Julia 在該基準測試上花費了 196 微秒，因此編譯的 Python 快 40%。
向量化
讓我們來看另一個示例。老實地講，我不確定要度量什么，但這是 Julia 團隊使用的代碼。
def parse_int(t):
for i in range(1,t):
n = random.randint(0,2**32-1)
s = hex(n)
if s[-1]=='L':
s = s[0:-1]
m = int(s,16)
assert m == n
return n
實際上，Julia 團隊的代碼有一條額外的指令，用于在存在末尾的 ‘L’ 時刪除它。我的 Anaconda 安裝需要這一行，但我的 Python 3安裝不需要它，所以我刪除了它。最初的代碼是：
def parse_int(t):
for i in range(1,t):
n = random.randint(0,2**32-1)
s = hex(n)
if s[-1]=='L':
s = s[0:-1]
m = int(s,16)
assert m == n
return n
對修改后的代碼計時會得到：
100 loops, best of 3:3.33 ms per loop
Numba 似乎沒什么幫助。Cython 代碼運行速度快了約 5 倍：
1000 loops, best of 3:617 μs per loop
Cython 代碼運行速度快了約 5 倍，但這還不足以彌補與 Julia 的差距。
我對此基準測試感到迷惑不解，我剖析了最初的代碼。以下是結果：
Timer unit:1e-06 s
Total time:0.013807 s
File:
Function: parse_int at line 1
Line # Hits Time Per Hit % Time Line Contents
==============================================================
1 def parse_int():
2 1000 699 0.7 5.1 for i in range(1,1000):
3 999 9149 9.2 66.3 n = random.randint(0,2**32-1)
4 999 1024 1.0 7.4 s = hex(n)
5 999 863 0.9 6.3 if s[-1]=='L':
6 s = s[0:-1]
7 999 1334 1.3 9.7 m = int(s,16)
8 999 738 0.7 5.3 assert m == n
可以看到，大部分時間都花費在了生成隨機數上。我不確定這是不是該基準測試的意圖。
加速此測試的一種方式是使用 Numpy 將隨機數生成移到循環之外。我們一次性創建一個隨機數數組。
def parse_int_vec():
n = np.random.randint(0,2**32-1,1000)
for i in range(1,1000):
ni = n[i]
s = hex(ni)
m = int(s,16)
assert m == ni
對它計時會得到：
1000 loops, best of 3:848 μs per loop
還不錯，快了 4 倍，接近于 Cython 代碼的速度。
擁有數組后，通過循環它來一次向某個元素應用 hex() 和 int() 函數似乎很傻。好消息是，Numpy 提供了一種向數組應用函數的方法，而不必使用循環，該函數是numpy.vectorize() 函數。此函數接受一次處理一個對象的函數。它返回一個處理數組的新函數。
vhex = np.vectorize(hex)
vint = np.vectorize(int)
def parse_int_numpy():
n = np.random.randint(0,2**32-1,1000)
s = vhex(n)
m = vint(s,16)
np.all(m == n)
return s
此代碼運行速度更快了一點，幾乎像 Cython 代碼一樣快：
1000 loops, best of 3:703 μs per loop
我肯定 Python 專家能夠比我在這里做得更好，因為我不太熟悉 Python 解析，但這再一次表明避免 Python 循環是個不錯的想法。
結束語
上面介紹了如何加快 Julia 團隊所使用的 4 個示例的運行速度。還有 3 個例子：
- pisum 使用 Numba 的運行速度快 29 倍。
- randmatstat 使用 Numpy 可將速度提高 2 倍。
- randmatmul 很簡單，沒有工具可應用到它之上。
包含所有 7 個示例的完整代碼的 Notebook 可在此處獲得。
我們在一個表格中總結一下我們的結果。我們給出了在最初的 Python 代碼與優化的代碼之間實現的加速。我們還給出了對 Julia 團隊使用的每個基準測試示例使用的工具。
這個表格表明，在前 4 個示例中，優化的 Python 代碼比 Julia 更快，后 3 個示例更慢。請注意，為了公平起見，對于 Fibonacci，我使用了遞歸代碼。
我認為這些小型的基準測試沒有提供哪種語言最快的明確答案。舉例而言， randmatstat 示例處理 5×5 矩陣。使用 Numpy 數組處理它有點小題大做。應該使用更大的矩陣執行基準測試。
我相信，應該在更復雜的代碼上對語言執行基準測試。Python 與 Julia 比較–一個來自機器學習的示例中提供了一個不錯的示例。在該文章中，Julia 似乎優于 Cython。如果我有時間，我會使用 Numba試一下。
無論如何，可以說，在這個小型基準測試上，使用正確的工具時，Python 的性能與 Julia 的性能不相上下。相反地，我們也可以說，Julia 的性能與編譯后的
Python 不相上下。考慮到 Julia 不需要對代碼進行任何注釋或修改，所以這本身就很有趣。
補充說明
我們暫停一會兒。我們已經看到在 Python 代碼性能至關重要時，應該使用許多工具：
- 使用 line_profiler 執行剖析。
- 編寫更好的 Python 代碼來避免不必要的計算。
- 使用向量化的操作和通過 Numpy 來廣播。
- 使用 Cython 或 Numba 編譯。
使用這些工具來了解它們在哪些地方很有用。與此同時，請謹慎使用這些工具。分析您的代碼，以便可以將精力放在值得優化的地方。重寫代碼來讓它變得更快，有時會讓它難以理解或通用性降低。因此，僅在得到的加速物有所值時這么做。Donald Knuth 曾經恰如其分地提出了這條建議：
“我們在 97% 的時間應該忘記較小的效率：不成熟的優化是萬惡之源。”
但是請注意，Knuth 的引語并不意味著優化是不值得的，例如，請查看停止錯誤地引用 Donald Knuth 的話!和‘不成熟的優化是惡魔’的謊言。
Python 代碼可以優化，而且應該在有意義的時間和位置進行優化。
來自：http://developer.51cto.com/art/201705/540479.htm