why do numpy so efficient

Posted on

之前在毕设中使用了numpy来完成进化算法的优化，比起使用python自带的数据类型list，速度提高了将近10倍，但是为什么numpy如此地高效，我在网上的博客看到的理解主要有两点：优化的存储管理和使用了低级语言完成性能优化。通过学习论文The NumPy Array: A Structure for Efficient Numerical Computation来详细了解一下其中的原理吧^_^。

numpy.ndarray structure:A View on Memory

自己的标题感觉没有论文的标题更生动，所以就借用了。numpy.ndarray是numpy的基本数据类型，它使用如下的数据结构来管理内存：

(1)data pointer: 存储ndarray在内存块中的首地址

(2)data type description: 数据类型

(3)shape: ndarray的shape属性

(4)strides: 存储为了维度搜索所需要跳过的byte数。这个只看定义有点抽象，可以通过下面的例子来理解：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
>>> import numpy as np
>>> a=np.random.random((10,10))
>>> a.dtype 
dtype('float64') # a每个元素是float64型的，占用8byte
>>> a.strides # 如果从a[0][1]跳转到a[1][0]则需要跳过10*8byte，从a[0][0]->a[0][1]则需要8byte
(80, 8)
>>> b=np.zeros((10,10)).astype(np.int)
>>> c=np.zeros((10,10)).astype(np.byte)
>>> b.dtype
dtype('int64') # b是int64型的，和a的strides应该是一样的
>>> b.strides
(80, 8) # 和c.strides一样
>>> c.dtype
dtype('int8') # c是int8型的，所以同理推理可以得到c.strides=(10,1)
>>> c.strides
(10, 1)

(5)flags: 定义的一些因素(factors),也是通过下面的例子看一下吧：

1
2
3
4
5
6
7
8
9
10
11
# 环境承接上面
# factors的解释来自numpy官网
>>> c=np.zeros((10,10)).astype(np.byte)
>>> c.flags
  C_CONTIGUOUS : True # c在内存中是按照C语言数据格式存储的，即按照行在内存中逐个存储
  F_CONTIGUOUS : False # c没有遵守Fortran的存储风格，即将列依次在内存中存储
  OWNDATA : True # True表示c占有这块内存，False表示内存是从另外的对象处借来的
  WRITEABLE : True # 代表这块内存区域是否是可写的，将它设置成False就可以是memory块变为read-only
  ALIGNED : True # 数据和所以的元素都和根据硬件环境适当地对其
  WRITEBACKIFCOPY : False # 后两个都是控制拷贝得到的array在改变数据时候对base array的操作
  UPDATEIFCOPY : False # 还有一些factor没有显示，可以参考官网

stride的机制给array访问相同部分的内存块提供了一种有效的方法，通过下面的例子可以看到：

1
2
3
4
5
6
7
8
9
>>> a=np.arange(9).reshape(3,3)
>>> a.strides
(24, 8)
>>> b=a[::2,::2]
>>> b
array([[0, 2],
       [6, 8]])
>>> b.strides # 0到6跨越了2行，所以对应stride是24*2=48
(48, 16)

a和b访问的是同一块内存块，我们称b是a的一个视图(view)。但是事实上创建one view of array不止可以通过切片操作，比如reshape，transpose，这些操作其实都是改变了array对象的strides，内存中的数据并没有发生变化，没有数据复制或者生成在这个过程中，所以这些view操作相当高效。下面给出transpose和reshape的例子：

1
2
3
4
5
6
7
8
9
>>> a=np.arange(12)
>>> b=a.reshape(3,4)
>>> c=b.T
>>> a.strides
(8,)
>>> b.strides
(32, 8)
>>> c.strides
(8, 32)

其实每个numpy.ndarray就是a view of memory，内存在view数据操作中并没有发生任何变化。

Numerical Operations on Arrays

在任何脚本语言(这一理念也适用于像C语言这样的低级语言吗)中，不当的使用for循环都会降低算法的性能，特别是对于大型数据集的每个元素应用简单计算(就像我毕设一开始使用的deap库。。。)，而一种提高性能的理念，就是通过矢量化(vectorization)的操作。其实这个思想在上matlab编程课的时候都老师也讲过，尽量通过矩阵运算来解决问题，而不是依赖for循环来操作。

numpy的矢量化编程中最重要的一个操作就是广播(broadcast)，这个term对我并不陌生，事实上我经常在numpy的Error中见到它T_T：

1
ValueError: operands could not be broadcast together with shapes (3,4) (3,)

但是broadcast的思想才是numpy矢量化操作的灵魂，我觉得它的魅力和numpy的花样索引一样吸引人，那么不废话，我们先看看官网对broadcast的说明:

broadcast描述了numpy在对不同维度的array进行数学运算时的操作。在了解broadcast之前，我在进行数学操作时都会保证array的维度匹配，这是最正常但是同时也是最耗时的操作，因为为了维度一致进行reshape，transpose这样的view操作还好，但是有的时候会进行一些数据的复制，对于大型数据来说会大大降低性能。如何用好broadcast，以及broadcast带来的性能提升有多少，我们先看下面一个例子：

1
2
3
4
5
6
7
8
9
# 我比较懒，这个标题下的代码都是复制的官网
>>> from numpy import array
>>> a = array([1.0, 2.0, 3.0])
>>> b = array([2.0, 2.0, 2.0])
>>> a * b
array([ 2.,  4.,  6.])
>>> b = 2.0
>>> a * b
array([ 2.,  4.,  6.])

第二次ab就用到了broadcast，我们可以简单地理解broadcast操作其实就是一个“expand”的过程，它将标量形式的b拓展成了和a维度一致的array。但是这个expand只是概念上的而不是物理上的，numpy is smart* enough to use the original scalar value without actually making copies so that broadcasting operations are as memory and computationally efficient as possible. 而且因为操作了更少的内存单元，第二种操作理论上会提高计算性能，测试表明对于一个array with 1 million values，broadcast可以提高10%的性能(但是现在的计算机估计性能提升就没这么明显了，这也就是为什么我的stupid code也可以跑的虎虎生风)。

对broadcast有了简明的认识之后，我们再来看看broadcast的一般规律：

broadcast rule: In order to broadcast, the size of the trailing axes for both arrays in an operation must either be the same size or one of them must be one.

numpy在计算中都会检测是否满足上面的broadcast rule，如果不满足就会出现开头处的ValueError。这个规律也很好理解，只要我们计算的array满足对应维度要么相等，要么有一方为1就行了，下面给出几个awesome use of broadcast:

1
2
3
4
5
6
7
8
9
>>> a=np.random.random([8,1,6,1])
>>> b=np.random.random([  7,1,5])
>>> (a+b).shape
(8, 7, 6, 5) # amazing!
# 经常会在图像处理中这样用
>>> image=np.random.random([256,256,3])
>>> scale=np.random.random([        3])
>>> (image*scale).shape
(256, 256, 3)

其实论文里面也讲了规则，但是没有例子来的生动，论文给了下面几个例子我觉得还挺有意思的：

use numpy in function

1
2
3
4
5
6
>>> def f(x):
...     return x**2-3*x+4
... 
>>> x=np.arange(1e5)
>>> y=[f(i) for i in x]
>>> y=f(x)

这个体现了numpy的矢量化编程理念，当你用for的时候，想一想能不能改进为矢量化操作呢？

using broadcasting to create a gird

这个看论文吧，就是broadcast操作的一种用法，注意如果使用broadcast的时候，你可以认为它是一种“expand”操作，因为这样好理解，你只要知道背后是使用了数学运算进行了优化就行。

computer vision

这一部分主要讲了numpy.dot()的使用，还有numpy.newaxis这个小知识点。dot就是矩阵的点乘，线性代数里我们经常用，但是如果dot的对象超过2维，也就是对tensor进行dot运算，结果遵循怎样的数学计算？通过查阅官网对应api我找到如下公式:

说实话看不太懂，等用到的时候再回来填坑吧。

关于numpy.newaxis的用法很简单，用来在array中插入一个新的axis，主要用来配合broadcast使用，下面给几个简单的例子：

1
2
3
4
5
6
7
8
9
10
11
12
>>> import numpy as np
>>> a=np.arange(9)
>>> b=a[np.newaxis,:]
>>> b
array([[0, 1, 2, 3, 4, 5, 6, 7, 8]])
>>> b.shape
(1, 9)
>>> a.shape
(9,)
>>> b=a[np.newaxis,np.newaxis,:]
>>> b.shape
(1, 1, 9)

Sharing Data

本来以为下午就写好，结果这会都要去吃饭了。。。晚上还要看tensorflow的论文，估计也是看不完了，不得不说numpy的特性真的多到人眼花缭乱T_T。(PS:配aria2将近1d,结果发现是网站不稳定，多试试就好了，真是醉了)。

晚上aria2出问题了，原因参考前一篇博客，现在我心态爆炸无心学习，这部分内容有空再补上吧，我要看看tensorflow的论文转换一下心情。这部分内容讲的是numpy如何高效利用外存，最常见的也就是磁盘，来存储数据，这里numpy使用了memory mapping方法。to be continuing …