优化代码
翻译自：http://scipy-lectures.github.com/advanced/optimizing/index.html
作者:Gaël Varoquaux
License:Creative Commons Attribution 3.0 United States License (CC-by) http://creativecommons.org/licenses/by/3.0/us
过早的优化是罪恶的根源。
——Donald. Knuth
这个章节涉及使Python代码运行更快的策略。
先决条件
- line_profiler(http://packages.python.org/line_profiler/)
目录
- toc
{: toc}
优化工作流
- 让它工作：以简单清晰的方式书写代码。
- 让它可靠的动作：书写自动化的测试实例，确认你的算法是正确的。如果你中止它，测试将捕捉到中断。
- 优化代码：通过剖析(profile)简单的用例来发现瓶颈，并且加速这些瓶颈，找到更好的算法或实现。记住在剖析一个现实的实例和代码的简洁与执行速度之间权衡。对于有效率的工作，最好让剖析运行约十秒。
剖析(profile)Python代码
非测量，不优化
- 测量：剖析(profiling)， 计时(timing)
Timeit
在Ipython中，使用timeit(http://docs.python.org/library/timeit.html)来计算基本运算时间。
In [1]: import numpy as np
In [2]: a = np.arange(1000)
In [3]: %timeit a ** 2
100000 loops, best of 3: 3.55 us per loop
In [4]: %timeit a ** 2.1
10000 loops, best of 3: 105 us per loop
In [5]: %timeit a * a
100000 loops, best of 3: 3.5 us per loop
使用这个指引你的策略选择。
注意：对于长时间运行的调用，使用%time代替%timeit;虽然精确度降低但运行更快。
分析器(Profiler)
当你有个很大的程序去分析时会有用，例如以下文件：
1 | mport numpy as np |
在IPython中我们可以测量脚本运行时间：
In [6]: %run -t demo.py
IPython CPU timings (estimated):
User : 11.03 s.
System : 0.43 s.
Wall time: 13.12 s.
然后分析(profile)它：
In [7]: %run -p demo.py
1169 function calls in 10.765 seconds
Ordered by: internal time
ncalls tottime percall cumtime percall filename:lineno(function)
2 10.693 5.346 10.699 5.350 decomp_svd.py:14(svd)
144 0.040 0.000 0.040 0.000 {numpy.core.multiarray.dot}
1 0.015 0.015 0.015 0.015 {method 'random_sample' of 'mtrand.RandomState' objects}
20 0.005 0.000 0.007 0.000 function_base.py:526(asarray_chkfinite)
1 0.003 0.003 10.764 10.764 demo.py:1(<module>)
18 0.002 0.000 0.002 0.000 decomp.py:197(eigh)
17 0.001 0.000 0.001 0.000 fastica_.py:219(gprime)
40 0.001 0.000 0.001 0.000 {method 'any' of 'numpy.ndarray' objects}
17 0.001 0.000 0.001 0.000 fastica_.py:215(g)
1 0.001 0.001 0.008 0.008 fastica_.py:88(_ica_par)
1 0.001 0.001 10.764 10.764 {execfile}
52 0.000 0.000 0.001 0.000 twodim_base.py:220(diag)
18 0.000 0.000 0.003 0.000 fastica_.py:40(_sym_decorrelation)
18 0.000 0.000 0.000 0.000 {method 'mean' of 'numpy.ndarray' objects}
显然svd(在decomp.py中)是占用我们时间最多的东西，即瓶颈。我们得找到一种方法来让这一步运行更快，或者避免这一步(算法优化)。加速代码剩下的部分并无益处。
Line-profiler
分析器(profiler)很棒：它告诉我们那个函数费时最多，但并不是它在哪里调用。
对此，我们使用line_profiler：在原文件中，我们用@profile修饰一些我们想要检查的函数(不用导入它)：
1 | @profile |
然后我们使用kernprof.py程序，用-l和-v：
lyy@arch ~ % 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 6
Total time: 10.5932 s
Line # Hits Time Per Hit % Time Line Contents
==============================================================
6 @profile
7 def test():
8 1 11070 11070.0 0.1 data = np.random.random((5000, 100))
9 1 10530291 10530291.0 99.4 u, s, v = linalg.svd(data)
10 1 31026 31026.0 0.3 pca = np.dot(u[:10, :], data)
11 1 20766 20766.0 0.2 results = fastica(pca.T)
kernprof.py -l -v demo.py 12.57s user 0.25s system 99% cpu 12.891 total
SVD占用了大部分时间。我们需要优化这行。
让代码更快
一旦我们确认了瓶颈，我们需要让相应的代码运行更快。
算法优化
首先寻找算法优化：有没有运算更少或更好的方式？
对于更高层次地看待问题，充分理解算法后的数学很有帮助。然而，寻找简单的改变并不寻常，像移动计算^1和在循环外分配内存^2,会带来很大益处。
####SVD的示例
在以上每个例子中，SVD——奇异值分解——是占用时间最多的。确实，当输入矩阵大小为n时算法计算代价大约是$n^3$。
然而，这些例子中，我们都没有使用SVD的输出，但是仅仅使用它最开始返回的参数最初很少的几行。如果我们使用scipy的svd实现，我们可以获得一个不完整的SVD版本。注意这个在scipy中的线性代数实现比numpy中的更加丰富，应该优先使用。
In [4]: %timeit np.linalg.svd(data)
1 loops, best of 3: 10.8 s per loop
In [5]: from scipy import linalg
In [6]: %timeit linalg.svd(data)
1 loops, best of 3: 10.4 s per loop
In [7]: %timeit linalg.svd(data, full_matrices=False)
1 loops, best of 3: 278 ms per loop
In [8]: %timeit np.linalg.svd(data, full_matrices=False)
1 loops, best of 3: 276 ms per loop
真正的不完全SVD，例如仅仅计算前十个特征向量，可以用arpack^3计算，可以在scipy.sparse.linalg.eigsh获得。
计算线性代数
对于确定的算法，许多瓶颈将是线性代数计算。在本例中，使用正确的函数解决正确的问题是关键。例如，一个对称矩阵的本征值问题比一个普通矩阵更容易解决。同样的，很多时候，你可以避免反转矩阵，并且使用代价更小(并更数值稳定)的运算。
了解你的线性代数计算。当有疑问时，探索scipy.linalg，并且用%timeit 来对你的数据尝试不同的选择。
写更快的数值代码
一个完整的关于使用numpy的讨论可以在Advanced Numpy章节中找到，或者在van der Walt等人的文章The NumPy array: a structure for efficient numerical computation中。这里我们仅仅讨论加速代码运行速度常见的技巧。
向量化循环
找到技巧来避免使用numpy数组循环。对此，掩码(masks)数组和索引(indices)数组会更有用。
广播
使用广播(broadcasting)来在结合它们之前对数组尽可能少的运算。
在适当的位置运算
In [9]: a = np.zeros(1e7) In [11]: %timeit global a ; a *= 0 10 loops, best of 3: 29.1 ms per loop in [12]: %timeit global a ; a = 0*a 10 loops, best of 3: 54.3 ms per loop注意：我们需要个
global a让timeit工作，因为它被赋给a，因此将它视作一个局部变量。善待内存：使用视图(views)而非拷贝(copies)
拷贝一个大数组就像对它们进行简单的数值计算一样耗费资源
In [18]: a = np.zeros(1e7) In [19]: %timeit a.copy() 10 loops, best of 3: 69 ms per loop In [20]: %timeit a + 1 10 loops, best of 3: 56.2 ms per loop小心缓存影响(cache effects)
内存存取当是成组时是省资源的：以连续的方式存取一个大数组比随机存取更快。这意味着除其它事项外更小的元素间距更快(参见CPU cache effects)^4
In [21]: c = np.zeros((1e4, 1e4), order=’C’)
In [22]: %timeit c.sum(axis=0) 1 loops, best of 3: 3.62 s per loop In [23]: %timeit c.sum(axis=1) 10 loops, best of 3: 171 ms per loop In [24]: c.strides Out[24]: (80000, 8) In [25]: c = np.zeros((1e4, 1e4), order='F') In [26]: %timeit c.sum(axis=0) 1 loops, best of 3: 166 ms per loop In [27]: %timeit c.sum(axis=1) 1 loops, best of 3: 3.63 s per loop这就是为何Fortran顺序或C顺序可能对运算的影响很大：
in [28]: a = np.random.rand(20, 2**18) in [29]: b = np.random.rand(20, 2**18) in [30]: %timeit np.dot(b, a.T) 1 loops, best of 3: 278 ms per loop in [31]: c = np.ascontiguousarray(a.T) in [32]: %timeit np.dot(b, c) 1 loops, best of 3: 1.94 s per loop
注意拷贝数据来解决这个影响不值得：
In [34]: %timeit c = np.ascontiguousarray(a.T)
10 loops, best of 3: 45.4 ms per loop
使用numexpr来自动为这种效应优化会很有用。
使用编译好的代码
一旦你确定所有高层次的优化都已经摸索过了，最后手段是将热点，也就是耗费时间最多的代码或函数，变成编译好的代码。对于编译代码，最优的选择是使用Cython:它很轻松地让你将已知的Python代码转换成编译好的代码，并且对numpy数组很好利用numpy支持产生有效率的代码，例如通过循环展开(unrolling loops)。
Waring:对上述所有流程，分析(profile)并且计时(time)你的选择。不要以理论为依据来进行优化。
{: alert }