2.4 代碼優化 · SciPy Lecture Notes 中文版

# 2.4 代碼優化 **作者**：Ga?l Varoquaux **Donald Knuth** "過早的優化是一切罪惡的根源" 本章處理用策略讓Python代碼跑得更快。 **先決條件** * line_profiler * gprof2dot * 來自dot實用程序 **章節內容** * 優化工作流 * 剖析Python代碼 * Timeit * Profiler * Line-profiler * Running `cProfile` * Using `gprof2dot` * 讓代碼更快 * 算法優化 * SVD的例子 * 寫更快的數值代碼 * 其他的鏈接 ## 2.4.1 優化工作流 1. 讓它工作起來：用簡單清晰的方式來寫代碼。 2. 讓它可靠的工作：寫自動的測試案例，以便真正確保你的算法是正確的，并且如果你破壞它，測試會捕捉到。 3. 通過剖析簡單的使用案例找到瓶頸，并且加速這些瓶頸，尋找更好的算法或實現方式來優化代碼。記住在剖析現實例子時簡單和代碼的執行速度需要進行一個權衡。要有效的運行，最好讓剖析工作持續10s左右。 ## 2.4.2剖析Python代碼 **無測量無優化！** * **測量**: 剖析, 計時 * 你可能會驚訝：最快的代碼并不是通常你想的樣子 ### 2.4.2.1 Timeit 在IPython中，使用timeit([http://docs.python.org/library/timeit.html)來計時基本的操作：](http://docs.python.org/library/timeit.html)來計時基本的操作：) In?[2]: ``` import numpy as np a = np.arange(1000) %timeit a ** 2 ``` ``` The slowest run took 60.37 times longer than the fastest. This could mean that an intermediate result is being cached 100000 loops, best of 3: 1.99 μs per loop ``` In?[3]: ``` %timeit a ** 2.1 ``` ``` 10000 loops, best of 3: 45.1 μs per loop ``` In?[4]: ``` %timeit a * a ``` ``` The slowest run took 12.79 times longer than the fastest. This could mean that an intermediate result is being cached 100000 loops, best of 3: 1.86 μs per loop ``` 用這個信息來指導在不同策略間進行選擇。 **筆記**：對于運行時間較長的單元，使用`%time`來代替`%timeit`; 它準確性較差但是更快。 ### 2.4.2.2 Profiler 當你有個大型程序要剖析時比較有用，例如[下面這個程序](http://scipy-lectures.github.io/_downloads/demo.py)： In?[?]: ``` # For this example to run, you also need the 'ica.py' file import numpy as np from scipy import linalg from ica import fastica def test(): data = np.random.random((5000, 100)) u, s, v = linalg.svd(data) pca = np.dot(u[:, :10].T, data) results = fastica(pca.T, whiten=False) if __name__ == '__main__': test() ``` **筆記**：這種技術是兩個非監督學習技術的組合，主成分分析（PCA）和獨立成分分析（[ICA](http://scipy-lectures.github.io/advanced/optimizing/index.html#id16)）。PCA是一種降維技術，即一種用更少的維度解釋數據中觀察到的變異的算法。ICA是一種源信號分離技術，例如分離由多個傳感器記錄的多種信號。如果傳感器比信號多，那么先進行PCA然后ICA會有幫助。更多的信息請見：[來自scikits-learn的FastICA例子](http://scikit-learn.org/stable/auto_examples/decomposition/plot_ica_blind_source_separation.html)。要運行它，你也需要下載[ica模塊](http://scipy-lectures.github.io/_downloads/ica.py)。在IPython我們計時這個腳本： In?[8]: ``` %run -t demo.py ``` ``` IPython CPU timings (estimated): User : 6.62 s. System : 0.17 s. Wall time: 3.72 s. ``` ``` /Users/cloga/Documents/scipy-lecture-notes_cn/ica.py:65: RuntimeWarning: invalid value encountered in sqrt W = (u * np.diag(1.0/np.sqrt(s)) * u.T) * W # W = (W * W.T) ^{-1/2} * W /Users/cloga/Documents/scipy-lecture-notes_cn/ica.py:90: RuntimeWarning: invalid value encountered in absolute lim = max(abs(abs(np.diag(np.dot(W1, W.T))) - 1)) ``` 并且剖析它： ``` %run -p demo.py 301 function calls in 3.746 seconds Ordered by: internal time ncalls tottime percall cumtime percall filename:lineno(function) 1 3.714 3.714 3.715 3.715 decomp_svd.py:15(svd) 1 0.019 0.019 3.745 3.745 demo.py:3(<module>) 1 0.007 0.007 0.007 0.007 {method 'random_sample' of 'mtrand.RandomState' objects} 14 0.003 0.000 0.003 0.000 {numpy.core._dotblas.dot} 1 0.001 0.001 0.001 0.001 function_base.py:550(asarray_chkfinite) 2 0.000 0.000 0.000 0.000 linalg.py:1116(eigh) 1 0.000 0.000 3.745 3.745 {execfile} 2 0.000 0.000 0.001 0.000 ica.py:58(_sym_decorrelation) 2 0.000 0.000 0.000 0.000 {method 'reduce' of 'numpy.ufunc' objects} 1 0.000 0.000 0.000 0.000 ica.py:195(gprime) 1 0.000 0.000 0.001 0.001 ica.py:69(_ica_par) 1 0.000 0.000 3.726 3.726 demo.py:9(test) 1 0.000 0.000 0.001 0.001 ica.py:97(fastica) 1 0.000 0.000 0.000 0.000 ica.py:192(g) 23 0.000 0.000 0.000 0.000 defmatrix.py:290(__array_finalize__) 4 0.000 0.000 0.000 0.000 twodim_base.py:242(diag) 1 0.000 0.000 3.746 3.746 interactiveshell.py:2616(safe_execfile) 10 0.000 0.000 0.000 0.000 {numpy.core.multiarray.array} 1 0.000 0.000 3.745 3.745 py3compat.py:279(execfile) 1 0.000 0.000 0.000 0.000 {method 'normal' of 'mtrand.RandomState' objects} 50 0.000 0.000 0.000 0.000 {isinstance} 10 0.000 0.000 0.000 0.000 defmatrix.py:66(asmatrix) 10 0.000 0.000 0.000 0.000 defmatrix.py:244(__new__) 9 0.000 0.000 0.000 0.000 numeric.py:394(asarray) 1 0.000 0.000 0.000 0.000 _methods.py:53(_mean) 1 0.000 0.000 0.000 0.000 {posix.getcwdu} 4 0.000 0.000 0.000 0.000 {method 'astype' of 'numpy.ndarray' objects} 6 0.000 0.000 0.000 0.000 defmatrix.py:338(__mul__) 2 0.000 0.000 0.000 0.000 linalg.py:139(_commonType) 4 0.000 0.000 0.000 0.000 {method 'view' of 'numpy.ndarray' objects} 1 0.000 0.000 0.000 0.000 posixpath.py:329(normpath) 5 0.000 0.000 0.000 0.000 {abs} 1 0.000 0.000 0.000 0.000 {open} 1 0.000 0.000 0.000 0.000 blas.py:172(find_best_blas_type) 1 0.000 0.000 0.000 0.000 blas.py:216(_get_funcs) 1 0.000 0.000 0.000 0.000 syspathcontext.py:64(__exit__) 3 0.000 0.000 0.000 0.000 {max} 6 0.000 0.000 0.000 0.000 {method 'transpose' of 'numpy.ndarray' objects} 1 0.000 0.000 0.000 0.000 posixpath.py:120(dirname) 2 0.000 0.000 0.000 0.000 linalg.py:101(get_linalg_error_extobj) 2 0.000 0.000 0.000 0.000 linalg.py:106(_makearray) 3 0.000 0.000 0.000 0.000 {numpy.core.multiarray.zeros} 6 0.000 0.000 0.000 0.000 defmatrix.py:928(getT) 1 0.000 0.000 0.000 0.000 syspathcontext.py:57(__enter__) 2 0.000 0.000 0.000 0.000 linalg.py:209(_assertNdSquareness) 7 0.000 0.000 0.000 0.000 {issubclass} 4 0.000 0.000 0.000 0.000 {getattr} 1 0.000 0.000 0.000 0.000 posixpath.py:358(abspath) 5 0.000 0.000 0.000 0.000 {method 'startswith' of 'unicode' objects} 2 0.000 0.000 0.000 0.000 linalg.py:198(_assertRankAtLeast2) 2 0.000 0.000 0.000 0.000 {method 'encode' of 'unicode' objects} 10 0.000 0.000 0.000 0.000 {method 'get' of 'dict' objects} 1 0.000 0.000 0.000 0.000 _methods.py:43(_count_reduce_items) 1 0.000 0.000 0.000 0.000 {method 'all' of 'numpy.ndarray' objects} 4 0.000 0.000 0.000 0.000 linalg.py:124(_realType) 1 0.000 0.000 0.000 0.000 syspathcontext.py:54(__init__) 1 0.000 0.000 0.000 0.000 posixpath.py:61(join) 1 0.000 0.000 3.746 3.746 <string>:1(<module>) 1 0.000 0.000 0.000 0.000 _methods.py:40(_all) 4 0.000 0.000 0.000 0.000 linalg.py:111(isComplexType) 2 0.000 0.000 0.000 0.000 {method '__array_prepare__' of 'numpy.ndarray' objects} 4 0.000 0.000 0.000 0.000 {min} 1 0.000 0.000 0.000 0.000 py3compat.py:19(encode) 1 0.000 0.000 0.000 0.000 defmatrix.py:872(getA) 2 0.000 0.000 0.000 0.000 numerictypes.py:949(_can_coerce_all) 6 0.000 0.000 0.000 0.000 {method 'append' of 'list' objects} 1 0.000 0.000 0.000 0.000 numerictypes.py:970(find_common_type) 1 0.000 0.000 0.000 0.000 {method 'mean' of 'numpy.ndarray' objects} 11 0.000 0.000 0.000 0.000 {len} 1 0.000 0.000 0.000 0.000 numeric.py:464(asanyarray) 1 0.000 0.000 0.000 0.000 {method '__array__' of 'numpy.ndarray' objects} 1 0.000 0.000 0.000 0.000 {method 'rfind' of 'unicode' objects} 2 0.000 0.000 0.000 0.000 {method 'upper' of 'str' objects} 1 0.000 0.000 0.000 0.000 posixpath.py:251(expanduser) 3 0.000 0.000 0.000 0.000 {method 'setdefault' of 'dict' objects} 1 0.000 0.000 0.000 0.000 {method 'diagonal' of 'numpy.ndarray' objects} 1 0.000 0.000 0.000 0.000 lapack.py:239(get_lapack_funcs) 1 0.000 0.000 0.000 0.000 {method 'rstrip' of 'unicode' objects} 1 0.000 0.000 0.000 0.000 py3compat.py:29(cast_bytes) 1 0.000 0.000 0.000 0.000 posixpath.py:52(isabs) 1 0.000 0.000 0.000 0.000 {method 'split' of 'unicode' objects} 1 0.000 0.000 0.000 0.000 {method 'endswith' of 'unicode' objects} 1 0.000 0.000 0.000 0.000 {sys.getdefaultencoding} 1 0.000 0.000 0.000 0.000 {method 'insert' of 'list' objects} 1 0.000 0.000 0.000 0.000 {method 'remove' of 'list' objects} 1 0.000 0.000 0.000 0.000 {method 'join' of 'unicode' objects} 1 0.000 0.000 0.000 0.000 {method 'index' of 'list' objects} 1 0.000 0.000 0.000 0.000 misc.py:126(_datacopied) 1 0.000 0.000 0.000 0.000 {sys.getfilesystemencoding} 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} ``` 很明顯`svd`（**decomp.py**中）占用了最多的時間，換句話說，是瓶頸。我們要找到方法讓這個步驟跑的更快，或者避免這個步驟（算法優化）。在其他部分花費時間是沒用的。 ### 2.4.2.3 Line-profiler profiler很棒：它告訴我們哪個函數花費了最多的時間，但是，不是它在哪里被調用。關于這一點，我們使用[line_profiler](http://packages.python.org/line_profiler/)：在源文件中，[[email?protected]](/cdn-cgi/l/email-protection)（不需要導入它）修飾了一些想要用檢查的函數： In?[?]: ``` @profile def test(): data = np.random.random((5000, 100)) u, s, v = linalg.svd(data) pca = np.dot(u[: , :10], data) results = fastica(pca.T, whiten=False) ``` 接著我們用[kernprof.py](http://packages.python.org/line_profiler/kernprof.py)來運行這個腳本，開啟`-l, --line-by-line`和`-v, --view`來使用逐行profiler，并且查看結果并保存他們： ``` kernprof.py -l -v demo.py Wrote profile results to demo.py.lprof Timer unit: 1e-06 s File: demo.py Function: test at line 5 Total time: 14.2793 s Line # Hits Time Per Hit % Time Line Contents ============================================================== 5 @profile 6 def test(): 7 1 19015 19015.0 0.1 data = np.random.random((5000, 100)) 8 1 14242163 14242163.0 99.7 u, s, v = linalg.svd(data) 9 1 10282 10282.0 0.1 pca = np.dot(u[:10, :], data) 10 1 7799 7799.0 0.1 results = fastica(pca.T, whiten=False) ``` SVD占用了幾乎所有時間，我們需要優化這一行。 ### 2.4.2.4 運行`cProfile` 在上面的IPython例子中，Ipython只是調用了內置的[Python剖析器](http://docs.python.org/2/library/profile.html)`cProfile`和`profile`。如果你想要用一個可視化工具來處理剖析器的結果，這會有幫助。 ``` python -m cProfile -o demo.prof demo.py ``` 使用`-o`開關將輸入剖析器結果到文件`demo.prof`。 ### 2.4.2.5 使用`gprof2dot` 如果你想要更加視覺化的剖析器輸入結果，你可以使用[gprof2dot](http://code.google.com/p/jrfonseca/wiki/Gprof2Dot)工具： In?[?]: ``` gprof2dot -f pstats demo.prof | dot -Tpng -o demo-prof.png ``` 這會生成下面的圖片： ![](http://scipy-lectures.github.io/_images/demo-prof.png) 這種方法打印了一個類似前一種方法的圖片。 ## 2.4.3 讓代碼更快一旦我們識別出瓶頸，我們需要讓相關的代碼跑得更快。 ### 2.4.3.1 算法優化第一件要看的事情是算法優化：有沒有計算量更小的方法或者更好的方法？從更高的視角來看這個問題，對算法背后的數學有一個很好的理解會有幫助。但是，尋找到像**將計算或內存分配移到循環外**這樣的簡單改變，來帶來巨大的收益，通常很困難。 #### 2.4.3.1.1 SVD的例子在上面的兩個例子中，SVD - [奇異值分解](http://en.wikipedia.org/wiki/Singular_value_decomposition) - 花費了最多的時間。確實，這個算法的計算成本大概是輸入矩陣大小的$n^3$。但是，在這些例子中，我們并不是使用SVD的所有輸出，而只是它第一個返回參數的前幾行。如果我們使用scipy的`svd`實現，我們可以請求一個這個SVD的不完整版本。注意scipy中的線性代數實現比在numpy中更豐富，應該被優選選用。 In?[20]: ``` %timeit np.linalg.svd(data) ``` ``` 1 loops, best of 3: 4.12 s per loop ``` In?[21]: ``` from scipy import linalg %timeit linalg.svd(data) ``` ``` 1 loops, best of 3: 3.65 s per loop ``` In?[22]: ``` %timeit linalg.svd(data, full_matrices=False) ``` ``` 10 loops, best of 3: 70.5 ms per loop ``` In?[23]: ``` %timeit np.linalg.svd(data, full_matrices=False) ``` ``` 10 loops, best of 3: 70.3 ms per loop ``` 接下來我們可以用這個發現來[優化前面的代碼](http://scipy-lectures.github.io/_downloads/demo_opt.py)： In?[24]: ``` import demo ``` In?[27]: ``` %timeit demo.test() ``` ``` 1 loops, best of 3: 3.65 s per loop ``` In?[28]: ``` import demo_opt ``` In?[29]: ``` %timeit demo_opt.test() ``` ``` 10 loops, best of 3: 81.9 ms per loop ``` 真實的非完整版SVD，即只計算前十個特征向量，可以用arpack來計算，可以在`scipy.sparse.linalg.eigsh`找到。 **計算線性代數** 對于特定的算法，許多瓶頸會是線性代數計算。在這種情況下，使用正確的方法來解決正確的問題是關鍵。例如一個對稱矩陣中的特征向量問題比通用矩陣中更好解決。同樣，更普遍的是，你可以避免矩陣逆轉，使用一些成本更低（在數字上更可靠）的操作。了解你的計算線性代數。當有疑問時，查找`scipy.linalg`，并且用`%timeit`來試一下替代方案。 ## 2.4.4 寫更快的數值代碼關于numpy的高級使用的討論可以在[高級numpy](http://scipy-lectures.github.io/advanced/advanced_numpy/index.html#advanced-numpy)那章，或者由van der Walt等所寫的文章[NumPy數組: 一種高效數值計算結構](http://hal.inria.fr/inria-00564007/en)。這里只是一些經常會遇到的讓代碼更快的小技巧。 * 循環向量化找到一些技巧來用numpy數組避免循環。對于這一點，掩蔽和索引通常很有用。 * 廣播在數組合并前，在盡可能小的數組上使用廣播。 * 原地操作 In?[30]: ``` a = np.zeros(1e7) %timeit global a ; a = 0*a ``` ``` 10 loops, best of 3: 33.5 ms per loop ``` ``` /Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/ipykernel/__main__.py:1: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future if __name__ == '__main__': ``` In?[31]: ``` %timeit global a ; a *= 0 ``` ``` 100 loops, best of 3: 8.98 ms per loop ``` **注意**: 我們需要在timeit中`global a`，以便正常工作，因為，向a賦值，會被認為是一個本地變量。 * 對內存好一點：使用視圖而不是副本復制一個大數組和在上面進行簡單的數值運算一樣代價昂貴： In?[32]: ``` a = np.zeros(1e7) ``` ``` /Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/ipykernel/__main__.py:1: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future if __name__ == '__main__': ``` In?[33]: ``` %timeit a.copy() ``` ``` 10 loops, best of 3: 28.2 ms per loop ``` In?[34]: ``` %timeit a + 1 ``` ``` 10 loops, best of 3: 33.4 ms per loop ``` * 注意緩存作用分組后內存訪問代價很低：用連續的方式訪問一個大數組比隨機訪問快很多。這意味著在其他方式中小步幅會更快（見[CPU緩存作用](http://scipy-lectures.github.io/advanced/advanced_numpy/index.html#cache-effects)）: In?[35]: ``` c = np.zeros((1e4, 1e4), order='C') ``` ``` /Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/ipykernel/__main__.py:1: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future if __name__ == '__main__': ``` In?[36]: ``` %timeit c.sum(axis=0) ``` ``` The slowest run took 5.66 times longer than the fastest. This could mean that an intermediate result is being cached 1 loops, best of 3: 80.9 ms per loop ``` In?[37]: ``` %timeit c.sum(axis=1) ``` ``` 10 loops, best of 3: 79.7 ms per loop ``` In?[38]: ``` c.strides ``` Out[38]: ``` (80000, 8) ``` 這就是為什么Fortran順序或者C順序會在操作上有很大的不同： In?[39]: ``` a = np.random.rand(20, 2**18) ``` In?[40]: ``` b = np.random.rand(20, 2**18) ``` In?[41]: ``` %timeit np.dot(b, a.T) ``` ``` 10 loops, best of 3: 23.8 ms per loop ``` In?[42]: ``` c = np.ascontiguousarray(a.T) ``` In?[43]: ``` %timeit np.dot(b, c) ``` ``` 10 loops, best of 3: 22.2 ms per loop ``` 注意，通過復制數據來繞過這個效果是不值得的： In?[44]: ``` %timeit c = np.ascontiguousarray(a.T) ``` ``` 10 loops, best of 3: 42.2 ms per loop ``` 使用[numexpr](http://code.google.com/p/numexpr/)可以幫助自動優化代碼的這種效果。 * 使用編譯的代碼一旦你確定所有的高級優化都試過了，那么最后一招就是轉移熱點，即將最花費時間的幾行或函數編譯代碼。要編譯代碼，優先選項是用使用[Cython](http://www.cython.org/)：它可以簡單的將Python代碼轉化為編譯代碼，并且很好的使用numpy支持來以numpy數據產出高效代碼，例如通過展開循環。 **警告**：對于以上的技巧，剖析并計時你的選擇。不要基于理論思考來優化。 ### 2.4.4.1 其他的鏈接 * 如果你需要剖析內存使用，你要應該試試[memory_profiler](http://pypi.python.org/pypi/memory_profiler) * 如果你需要剖析C擴展程序，你應該用[yep](http://pypi.python.org/pypi/yep)從Python中試著使用一下[gperftools](http://code.google.com/p/gperftools/?redir=1)。 * 如果你想要持續跟蹤代碼的效率，比如隨著你不斷向代碼庫提交，你應該試一下：[vbench](https://github.com/pydata/vbench) * 如果你需要一些交互的可視化為什么不試一下[RunSnakeRun](http://www.vrplumber.com/programming/runsnakerun/)