@Copyright悟空实验室spacegoing2014

转载请注明作者、出处

High Performance Python

本文参考《Python High Performance

Programming》，简述如何编写科学计算领域高性能Python代码，涉及ProfileNumpy Numexpr

Cython和Library 中OptionalOperating System

Services部分内容。对于科学计算，强烈推荐使用Anaconda发行版，Spyder作为IDE并使用Ipythonconsole。本文对ipdb、conda、pypi的使用有简要描述。

友情提示：鉴于目前中文资料、社区太烂，这个层次面对的问题如果出错往往只有dev

team能够解决，英文稍逊的读者不建议使用Python进行高性能计算。

本文非入门文章，适合有Python基础和Matlab使用经验、追求性能的开发人员阅读。

除第四部分，本文以2.7版本为讲解，代码可运行于python2.7和python

3.4。第四部分仅可在3.4运行，有需要的读者可以参考library进行简要修改。

将近两周的时间被Python虐的很惨，为平息怒火使文章尽量中肯，首先简短谈下使用Python两周以来感受到的优点：灵活。抛去陈词滥调的Module、DataStructure等话题，开源是Python灵活的源泉，你可以随意修改底层代码定制属于你的Python，想象一下给MatlabPCT的parfor结构添加Initializer。配合强大的社区，这种灵活让你可以实现任何功能。

除此之外，与Matlab相比，对科学计算任务而言，Python弱爆了。(网上充斥的观点是Matlab有多慢多慢，而长时间以来我亲身实践表明说这些话的人根本不会用Matlab。你真的了解矩阵运算吗？你认真看过你使用的工具箱的ReleaseNote各个版本间的优化说明吗？如果你还在用r2014a以下的版本，别看这篇文章了)。事实是，Matlab已经对我们这个层次能用到的所有运算进行了加速，除非你是大牛，否则自己用C写的算法往往不如Matlab快。

废话不多说，如果碰巧你需要的功能Matlab无法实现的话，按照以下步骤能让你科学计算代码的运算速度更加贴近Matlab。

(一) Ipython:Profiling & Ipdb

( Ipython Debug )

充分了解你代码的每个细节是提高代码性能的前提。优化代码是一个逐层细化的过程：首先找到需要优化的函数，然后找到某个函数中需要优化的某个方法，再明确某个方法中需要修改的代码行。

直接告诉大家我目前找到的最方便的方法：Ipython。Ipython集成了各种方便开发的Magics(Ipython中的jargon，其实就是一些小工具)，这里主要讲述%timeit、%prun、%lprun。由浅入深讲述如何profile整段代码、整段代码中的各个函数、以及函数中的每个表达式。

1.1 %timeit(相当于matlab中tic toc，但不能类似t=toc赋值)

%timeit Magic其实源于Python的timeit

Module。Timeit可以测试运行一段代码片段(snippets)消耗的时间，有两个输入参数：n指定每次测试运行多少次；r指定一共进行多少次测试。Ipython的%timeit可以根据snippets的计算量自动调整n和r，默认值n=100r=3。

比如有如下代码，分别用三种方法，对[1000000,1]的矩阵赋值：

#%%

data=[]

for i in range(n):

data.append(i)

def testapp(n=int(1e6)):

data=[]

apd=data.append

for i in range(n):

apd(i)

data=[None]*n

for i in range(n):

data[i]=i

#%%

def test():

testap()

testapp()

testappp()

在Ipython console中执行后(或在Spyder中按F5)，键入代码如下：

In [3]: %timeit test()

Out [3]: 1 loops, best of 3: 221ms per loop

可知，运行test()消耗的时间是221ms。

小知识1：代码段中的#%%为Spyder中特有的Token，将整个代码分成了两段cells(相当于Matlab中的%%把代码分成Sections)，使用ctrl+Enter可以单独运行当前选中的cell，使用ctrl+shift+enter可以在运行完毕当前cell后转入下一cell。

1.2 %prun: Profiling each call

很多情况下我们不会满足于测试某函数本身，更多的我们想知道函数内部，调用的所有函数的，时间消耗情况，这时我们就可以使用%prunMagic：

In [4]: %prun test()

Out [4]: Ordered by: internaltime

ncalls tottime percall cumtime percallfilename:lineno(function)

1 0.141 0.141 0.209 0.209

:1(testap)

2000000 0.126 0.000 0.126 0.000 {method 'append' of'list' objects}

1 0.078 0.078 0.135 0.135

:7(testapp)

1 0.051 0.051 0.061 0.061

:14(testappp)

1 0.020 0.020 0.425 0.425

:1(test)

1 0.010 0.010 0.010 0.010

{range}

1 0.000 0.000 0.425 0.425

:1()

1 0.000 0.000 0.000 0.000

{method 'disable' of '_lsprof.Profiler' objects}

这时我们就可以很轻松的看到函数内部的调用情况，每个指标含义如下：

Ncalls：函数调用次数

Tottime/percall：单独调用这个函数所消耗的总时间，以及每次调用平均消耗的时间

Cumtime/percall：调用这个函数(包括这个函数内调用的其它函数)消耗掉的时间，以及每次调用平均消耗的时间。

Filenames:lineno(function) ：文件名：行数(函数名)

小知识2：Pre-allocating

通过观察testap

testapp 和 testappp的执行时间我们可以发现Python中仍然需要Pre-allocating。官方文档(PythonReference)中写到每次append的计算复杂度是O(1)，但没有包括调用方法消耗的时间。Testapp中将Data.append方法赋值给变量apd，直接执行apd的时间比testap中的时间显著降低。

但即使如此，Pre-allocating的testappp显然还是王道。

1.3 %lprun： Profiling Line-by-Line

High

Performance来源于追求极致，来源于对每个变量的斤斤计较，显然对每个函数的profile仍然不能满足我们的需求，%lprunMagic就是为了那些逐句优化的苦逼设计的：

In [5]: %lprun -f testappp test()

Out [5]: Timer unit: 2.8428e-07 s

Total time: 0.789642 s

File:

Function: testappp at line 14

Line

# Hits Time Per Hit %

Time Line Contents

==============================================================

14 deftestappp(n=int(1e6)):

15 1 12752 12752.0 0.5 data=[None]*n

16 1000001 1247773 1.2 44.9 for i in range(n):

17 1000000 1517163 1.5 54.6 data[i]=i

每个指标含义如下：

Line#: 行数

Hits: 调用次数

Time per Hit: 每次调用消耗的时间

%Time: 占总时间的百分比

我们看到，虽然经过Pre-allocating的testappp()已经最快，但其中forin

range(n)居然占用了总执行时间的45%，这显然不正常。原因是与3.4不同，在Python2.7中，range(n)返回一个列表而非generator。因此在2.7中所有的循环都应该用xrange(返回generator)实现，将代码修改如下：

def testappp(n=int(1e6)):

data=[None]*n

for i in xrange(n):

data[i]=i

重新执行%lprun我们得到如下结果：

Total time: 0.688533 s

File:

Function: testappp at line 14

Line

# Hits Time Per Hit %

Time Line Contents

==============================================================

14 deftestappp(n=int(1e6)):

15 1 12798 12798.0 0.5 data=[None]*n

16 1000001 1105483 1.1 35.6 for i in xrange(n):

17 1000000 1303742 1.3 63.8 data[i]=i

我们看到无论总时间还是占比都显著降低了。

小知识3：使用Conda + Python PackagesIndex(PyPi)安装Module

有些读者可能在执行%lprun时会报错，有两种情况，如果你使用Python3.4，line_profiler不支持，你可以选择使用Spyder自带的line_profiler；如果报错找不到Module，那么参考如下指南。

Python灵活性除语言本身外，很大的原因就是社区支持。PyPi集成了网络上能够找到的绝大部分Module资源。但纷繁发杂的Module及其安装方法让人眼花缭乱。令人欣慰的是annaconda推出的conda架构大大简化了这个过程。这里仅以安装line_profiler为例，需要安装line_profiler，请打开AnacondaCommand

Prompt后输入如下命令：

pip install line_profiler

安装即可完成。

(运行前确保当前生效的是python2.7运行环境，详细内容请参考第四部分对Conda的讲解)

安装完成后，现在Ipython内加载扩展Module：

In [1]: %load_ext line_profiler

此时即可使用%lprun测试代码。

1.4 Ipython Debug

特别值得一提的是Ipython Debug。ipdb与gdb操作相似，但对熟悉GUIF5 F9

debug的Matlab用户这里还是再说明一下。

比如运行如下代码：

def testap(n=int(1e6)):

#data=[]

for i in range(n):

data.append(i)

test()

Ipython解释器会报错:

Traceback (most recent call last):

File "", line 1, in

test()

File "", line 2, in test

testap()

File "", line 5, in testap

data.append(i)

NameError: global name 'data' is notdefined

在Ipython

Console键入命令 �bug 直接进入ipdb调试模式：

In [2]: �bug

Out [2]: > (5)testap()

3

# data=[]

4 for i in xrange(n):

---->

5 data.append(i)

ipdb>

在调试模式中可以使用如下7个命令：

up: 返回调用本层函数的上层函数

down: 返回本层函数调用的下层函数

s: 执行选中行，如果使用了函数则ipdb进入函数体内继续debug

n: 执行选中行，如果使用了函数直接调用函数并得到返回值，不进入函数

r: 执行当前函数中的全部代码，在执行当前函数的return命令前停止执行;

或者在下一断点处停止

p: 打印p后的变量，比如 p

i 则在console中打印变量i。

exit: 退出debug模式

在ipdb中你仍然可以对变量赋值on the

flow，ipdb功能还是比较全面的，结合Spyder对于Matlab用户上手应该不太困难。

(二) Numpy与Numexpr

2.1 Numpy Broadcast简易教程

Python原生数据类型并不是为科学计算特殊设计的。笨拙的list和dict提高通用性同时不免以牺牲性能为代价。Scipy.org开发的Numpy项目专门为科学计算提供了解决方案。

限于时间精力有限，关于Numpy的详细教程请参考官方document，对于Matlab用户，特别推荐scipy.org编写的说明：http://wiki.scipy.org/NumPy_for_Matlab_Users，不夸张的说熟悉Matlab的用户可以很轻松地在1小时内学会Numpy的基本操作。

为了能让你在30分钟内学会(结合上面的连接)，这里特别讲一下Numpy与Matlab最大的不同：Broadcast机制。

Matlab用户不会对这个陌生，所谓的Broadcast其实就相当于Matlab中的scalar expansion，即：

A=zeros(5,1);

A=20;

这时你会发现5行1列的数组A中的每个元素都是20。

Broadcast正是这个概念：如果操作符(如”=”)两边大小不一致，那么，numpy将会自动选取变量的某一axis(jargon

in

Numpy，指变量的维数)扩展小的那个，并进行element-wise(针对每个元素)的操作。这里以两个例子进行说明。

首先看下面这个例子：

Import numpy as np

X=np.random.rand(5,1)

Y=5

D=X-Y

这个例子的结果是对X的每个元素都减去了5。原因很好理解，X的大小为[5,1]，Y的大小为[1,1]。因此在进行X-Y操作时，Numpy选取axis=0也就是X的第一维Broadcast

Y，就相当于Y是一个[5,1]的所有元素都为5的数组，然后再进行element-wise操作，即得我们想要的结果。

那么我们挑战一下极限，如果操作符两边变量大小都不为1，且不一致，这时Numpy会怎么办呢？比如经典的欧式距离矩阵问题，设现在二维空间中有5个粒子，并且知道这些粒子的坐标，现在要求5个粒子间两两的距离。首先随机生成5个粒子:

import numpy as np

dist=np.empty((n,n))

m=0

for i in particles:

n=0

for j in particles:

x_ij=i[0]-j[0]

y_ij=i[1]-j[1]

dist[m,n]=np.sqrt(x_ij**2+y_ij**2)

n+=1

m+=1

很自然地想到可以用两层、25次循环解决问题，生成一个[5,5]的距离矩阵保存这些数据。但是等等，生成[5,5]的距离矩阵？？！！

之前我们已经说过，broadcast的作用是将扩展一个变量使其和另外一个变量的大小一致，然后进行element-wise的操作。而距离矩阵的本质是在横轴方向和纵轴方向，分别扩展了5个粒子，先求出他们在xy方向上的绝对距离，然后开方求解欧式距离问题。那么我们是否可以利用Broadcast帮助我们减少某些循环，或者减少某些循环中的操作呢？

设想我们有粒子1到粒子5，我们首先将他们分别按照横、纵轴方向排列如下，然后对每个元素，两两相减、并求欧式距离，生成一个[5,5,2]的矩阵即可：

表2.1 粒子距离矩阵

Particle #1

Particle #2

Particle #3

Particle #4

Particle #5

Particle #1

[X1 Y1] [X1 Y1]

[X2 Y2] [X1 Y1]

[X3 Y3] [X1 Y1]

[X4 Y4] [X1 Y1]

[X5 Y5] [X1 Y1]

Particle #2

[X1 Y1] [X2 Y2]

[X2 Y2] [X2 Y2]

[X3 Y3] [X2 Y2]

[X4 Y4] [X2 Y2]

[X5 Y5] [X2 Y2]

Particle #3

[X1 Y1] [X3 Y3]

[X2 Y2] [X3 Y3]

[X3 Y3] [X3 Y3]

[X4 Y4] [X3 Y3]

[X5 Y5] [X3 Y3]

Particle #4

[X1 Y1] [X4 Y4]

[X2 Y2] [X4 Y4]

[X3 Y3] [X4 Y4]

[X4 Y4] [X4 Y4]

[X5 Y5] [X4 Y4]

Particle #5

[X1 Y1] [X5 Y5]

[X2 Y2] [X5 Y5]

[X3 Y3] [X5 Y5]

[X4 Y4] [X5 Y5]

[X5 Y5] [X5 Y5]

表2.2 粒子X\Y轴绝对距离矩阵

Particle #1

Particle #2

Particle #3

Particle #4

Particle #5

Particle #1

[X1-X1] [Y1-Y1]

[X2-X1] [Y2-Y1]

[X3-X1] [Y3-Y1]

[X4-X1] [Y4-Y1]

[X5-X1] [Y5-Y1]

Particle #2

[X1-X2] [Y1-Y2]

[X2-X2] [Y2-Y2]

[X3-X2] [Y3-Y2]

[X4-X2] [Y4-Y2]

[X5-X2] [Y5-Y2]

Particle #3

[X1-X3] [Y1-Y3]

[X2-X3] [Y2-Y3]

[X3-X3] [Y3-Y3]

[X4-X3] [Y4-Y3]

[X5-X3] [Y5-Y3]

Particle #4

[X1-X4] [Y1-Y4]

[X2-X4] [Y2-Y4]

[X3-X4] [Y3-Y4]

[X4-X4] [Y4-Y4]

[X5-X4] [Y5-Y4]

Particle #5

[X1-X5] [Y1-Y5]

[X2-X5] [Y2-Y5]

[X3-X5] [Y3-Y5]

[X4-X5] [Y4-Y5]

[X5-X5] [Y5-Y5]

那么问题来了，可以用挖掘机控制键盘写程序吗？我们注意到Particles是一个[2,5]的矩阵。文章开头就提到过，Broadcast的含义是，numpy将会自动选取变量的某一axis(jargonin

Numpy，指变量的维数)扩展小的那个，并进行element-wise(针对每个元素)的操作。那么如何将Particles扩展为一个[5,5,2]的矩阵呢？这种情况貌似没有任何一维可以Broadcast？

我们可以注意到，表2.2的[5,5,2]矩阵实际上是由Particles矩阵，在横、纵轴两个方向上扩展而来的(表2.1对应相减即为表2.2)。那么我们接下来需要做的事情就非常简单：构造一个沿X轴分布的Particles数组，再构建一个沿Y轴分布的数组即可。

X：

Particle #1

Particle #2

Particle #3

Particle #4

Particle #5

[X1 Y1]

[X2 Y2]

[X3 Y3]

[X4 Y4]

[X5 Y5]

Y:

Particle #1

[X1 Y1]

Particle #2

[X2 Y2]

Particle #3

[X3 Y3]

Particle #4

[X4 Y4]

Particle #5

[X5 Y5]

聪明如你们肯定想到了，我们期望生成的表2.2是一个[5,5,2]的矩阵，那么我们只需要构造一个[1,5,2]的X矩阵，在构造一个[5,1,2]的Y矩阵，那么Numpy就可以将X沿纵轴方向扩展、将Y沿横轴方向扩展了。

为了解决这个问题，numpy专门为我们提供了一个工具：np.newaxis，代码如下：

import numpy as np

def dist_matrix(n=5):

particles=np.random.rand(n,2)

X=particles[:,np.newaxis,:]

Y=particles[np.newaxis,:,:]

adist=X-Y

odist=np.sqrt((adist**2).sum(axis=2))

我们比较一下两个函数的执行时间：

In [2]: %timeit dist_loop()

Out [2]: 10000 loops, best of 3: 170µs per loop

In [3]: %timeit dist_matrix()

Out [3]: 10000 loops, best of 3:19.2 µs per loop

一定不需要我再多说什么了。

至此Numpy就已经讲述的差不多了，读者需要在了解一下numpy的Indexing机制(称作fancyIndex)，并对照Numpyfor

Matlab

User’s，了解Numpy各个数据类型的使用机制。余下的工作就是参考Scipy.org官方的docs对Numpy的各个工具进行了解。读者在需要特定功能时查阅即可。

小知识4:

Numpy的axis是从0开始计算的，我们说的第一维实际上是axis=0。Dist_matrix代码odist=np.sqrt((adist**2).sum(axis=2))中的axis=2意味着沿第三维求和。Adist的大小是[10,10,2]，沿第三维求和后odist的大小变为[10,10]正是我们要的结果。

小知识5：

特别需要注意np.empty的使用。此方法用于对变量进行preallocating。但是需要注意的是，empty并不是真正的None，此方法会给变量随机分配一块内存，变量中存储的值是之前残留的值，是完全随机的。因此很有可能分配的值是nan。因此如果你有类似操作:

Data1=np.empty((n,1))

Data2=np.empty((n+1,2))

for i inxrange(n):

Data2[i+1]=Data2[i]+Data1[i]

那么特别需要注意的是一定要先给Data2[0]赋值成一个对你的操作identity的值。比如此处的操作为”+”，那么Data2[0]=0。如果为’*’那么Data2[0]=1。否则你程序的结果将随机取决于分配的结果。

2.2 Numexpr

numexpr是由David M.

Cooke编写的，经过CPU缓存优化的Module。Numpy在进行复杂运算时会将中间结果储存在CPU缓存中，这不仅浪费缓存而且会增加操作从而减慢运算速度。numexpr正是从这方面着手解决了缓存优化问题。

使用numexpr非常简单，这里给出两点必须事项和一点建议：

必须：

1. 不能在numexpr中使用任何index和slice操作。

2. 如果表达式中具有reduction

operation(归约操作)，则必须将reductionoperation放置在最外层。

建议：

在一条语句中完成尽可能多的运算。

下面我们就是用numexpr来修改刚刚编写的numpy的代码dist_matrix()：

import numpy as np

import numexpr as ne

def dist_numexpr(n=1e4):

particles=np.random.rand(n,2)

X=particles[:,np.newaxis,:]

Y=particles[np.newaxis,:,:]

dist=ne.evaluate('sum((X-Y)**2,2)') #一个表达式中完成尽可能多的运算，但reduxtionoperation只能在最外层

dist=ne.evaluate('sqrt(dist)')

接下来我们看一下运行时间：

In [3]: %timeit dist_numexpr()

Out [3]: 1 loops, best of 3: 1.08s per loop

In [4]: %timeit dist_matrix(1e4)

Out [4]: 1 loops, best of 3: 4.13s per loop

尽管numexpr显著的降低了运行时间，但是在数据量较大时可能造成CPU缓存溢出的问题。