[转]numpy性能优化

转自：http://blog.csdn.net/pipisorry/article/details/39087583

http://blog.csdn.net/pipisorry/article/details/39087583

Introduction

NumPy提供了一个特殊的数据类型ndarray，其在向量计算上作了优化。这个对象是科学数值计算中大多数算法的核心。

相比于原生的Python，利用NumPy数组能够得到显著的性能加速，尤为是当你的计算遵循单指令多数据流（SIMD）范式时。

然而，利用NumPy也有可能有意无心地写出未优化的代码。下面这些技巧能够帮助你编写高效的NumPy代码。

避免没必要要的数据拷贝

查看数组的内存地址

1. 查看静默数组拷贝的第一步是在内存中找到数组的地址。下边的函数就是作这个的：

def  id (x):      # This function returns the memory block address of an array.      return x.__array_interface__[ 'data' ][ 0 ]

2. 有时你可能须要复制一个数组，例如你须要在操做一个数组时，内存中仍然保留其原始副本。

a  = np.zeros( 10 ); aid = id (a); aid 71211328 b  = a.copy(); id (b) = = aid False

Note:python列表的浅拷贝和深拷贝-列表拷贝一节

具备相同数据地址（好比id函数的返回值）的两个数组，共享底层数据缓冲区。然而，共享底层数据缓冲区的数组，只有当它们具备相同的偏移量（意味着它们的第一个元素相同）时，才具备相同的数据地址。共享数据缓冲区，但偏移量不一样的两个数组，在内存地址上有细微的差异：

id (a), id (a[ 1 :]) ( 71211328 , 71211336 )

在这篇文章中，咱们将确保函数用到的数组具备相同的偏移量。

下边是一个判断两个数组是否共享相同数据的更可靠的方案：

def get_data_base(arr):      "" "For a given Numpy array, finds the base array that " owns " the actual data." ""      base = arr      while is instance(base.base, np.ndarray):          base = base.base      return  base   def arrays_share_data(x, y):      return  get_data_base(x) is  get_data_base(y)   print(arrays_share_data(a,a.copy()), arrays_share_data(a,a[ 1 :])) False True

感谢Michael Droettboom指出这种更精确的方法，提出这个替代方案。

Note:a和a[1:]id虽然不一样，可是他们是共享内存的。

就地操做和隐式拷贝操做

3. 数组计算包括就地操做（下面第一个例子：数组修改）或隐式拷贝操做（第二个例子：建立一个新的数组）。

1 2 3 4 5

a  * = 2 ; id (a) = = aid True   c  = a * 2 ; id (c) = = aid False

必定要选择真正须要的操做类型。隐式拷贝操做很明显很慢，以下所示：

1 2 3 4 5 6 7

% % timeit a = np.zeros( 10000000 ) a  * = 2 10  loops, best of  3 : 19.2 ms per loop   % % timeit a = np.zeros( 10000000 ) b  = a * 2 10  loops, best of  3 : 42.6 ms per loop

4. 重塑数组可能涉及到拷贝操做，也可能涉及不到。

例如，重塑一个二维矩阵不涉及拷贝操做，除非它被转置（或更通常的非连续操做）：

1 2	a  = np.zeros(( 10 , 10 )); aid  =  id (a); aid 53423728

重塑一个数组，同时保留其顺序，并不触发拷贝操做。

1 2	b  = a.reshape(( 1 , - 1 )); id (b) = = aid True

转置一个数组会改变其顺序，因此这种重塑会触发拷贝操做。

1 2	c  = a.T.reshape(( 1 , - 1 )); id (c) = = aid False

所以，后边的指令比前边的指令明显要慢。

5. 数组的flatten和revel方法将数组变为一个一维向量（铺平数组）。flatten方法老是返回一个拷贝后的副本，而revel方法只有当有必要时才返回一个拷贝后的副本（因此该方法要快得多，尤为是在大数组上进行操做时）。

1 2 3 4 5 6 7 8 9 10 11

d  = a.flatten(); id (d) = = aid False   e  = a.ravel(); id (e) = = aid True   % timeit a.flatten() 1000000  loops, best of  3 : 881 ns per loop   % timeit a.ravel() 1000000 loops, best of 3 : 294 ns per loop

广播规则

广播规则容许你在形状不一样但却兼容的数组上进行计算。换句话说，你并不老是须要重塑或铺平数组，使它们的形状匹配。

广播规则描述了具备不一样维度和/或形状的数组仍能够用于计算。通常的规则是：当两个维度相等，或其中一个为1时，它们是兼容的。NumPy使用这个规则，从后边的维数开始，向前推导，来比较两个元素级数组的形状。最小的维度在内部被自动延伸，从而匹配其余维度，但此操做并不涉及任何内存复制。

下面的例子说明了两个向量之间进行矢量积的两个方法：第一个方法涉及到数组的变形操做，第二个方法涉及到广播规则。显然第二个方法是要快得多。

 

n =1000

a =np.arange(n)
ac =a[:, np.newaxis]
ar =a[np.newaxis, :]

%timeit np.tile(ac, (1, n))* np.tile(ar, (n,1))
100 loops, best of 3:10 ms per loop

%timeit ar* ac
100 loops, best of 3:2.36 ms per loop

[numpy教程- 通用函数ufunc- 广播规则]

 

NumPy数组进行高效的选择

NumPy提供了多种数组分片的方式。

数组视图涉及到一个数组的原始数据缓冲区，但具备不一样的偏移量，形状和步长。NumPy只容许等步长选择（即线性分隔索引）。

NumPy还提供沿一个轴进行任意选择的特定功能。

最后，花式索引（fancy indexing）是最通常的选择方法，但正如咱们将要在文章中看到的那样，它同时也是最慢的。

1. 建立一个具备不少行的数组。咱们将沿第一维选择该数组的分片。

n, d  = 100000 , 100 a  = np.random.random_sample((n, d)); aid = id (a)

数组视图和花式索引

2. 每10行选择一行，这里用到了两个不一样的方法（数组视图和花式索引）。

b1  = a[:: 10 ] b2  = a[np.arange( 0 , n, 10 )] np.array_equal(b1, b2) True

3. 数组视图指向原始数据缓冲区，而花式索引产生一个拷贝副本。

id (b1) = = aid, id (b2) = = aid ( True , False )

两个方法的执行效率，花式索引慢好几个数量级，由于它要复制一个大数组。

替代花式索引：索引列表

当须要沿一个维度进行非等步长选择时，数组视图就无能为力了。

然而，替代花式索引的方法在这种状况下依然存在。给定一个索引列表，NumPy的函数能够沿一个轴执行选择操做。

i  = np.arange( 0 , n, 10 )   b1  = a[i] b2  = np.take(a, i, axis = 0 )   np.array_equal(b1, b2) True

第二个方法更快一点：

1 2 3 4 5

% timeit a[i] 100  loops, best of  3 : 13 ms per loop   % timeit np.take(a, i, axis = 0 ) 100  loops, best of  3 : 4.87 ms per loop

替代花式索引：布尔掩码

当沿一个轴进行选择的索引是经过一个布尔掩码向量指定时，compress函数能够做为花式索引的替代方案。

1	i  = np.random.random_sample(n) < . 5

可使用花式索引或者np.compress函数进行选择。

b1  = a[i] b2  = np.compress(i, a, axis = 0 )   np.array_equal(b1, b2) True   % timeit a[i] 10  loops, best of  3 : 59.8 ms per loop   % timeit np.compress(i, a, axis = 0 ) 10  loops, best of  3 : 24.1 ms per loop

第二个方法一样比花式索引快得多。

花式索引是进行数组任意选择的最通常方法。然而，每每会存在更有效、更快的方法，应尽量首选那些方法。

当进行等步长选择时应该使用数组视图，但须要注意这样一个事实：视图涉及到原始数据缓冲区。

 

 

为何NumPy数组如此高效？

一个NumPy数组基本上是由元数据（维数、形状、数据类型等）和实际数据构成。数据存储在一个均匀连续的内存块中，该内存在系统内存（随机存取存储器，或RAM）的一个特定地址处，被称为数据缓冲区。这是和list等纯Python结构的主要区别，list的元素在系统内存中是分散存储的。这是使NumPy数组如此高效的决定性因素。

为何这会如此重要？主要缘由是：

1. 低级语言好比C，能够很高效的实现数组计算（NumPy的很大一部分其实是用C编写）。例如，知道了内存块地址和数据类型，数组计算只是简单遍历其中全部的元素。但在Python中使用list实现，会有很大的开销。

2. 内存访问模式中的空间位置访问会产生显著地性能提升，尤为要感谢CPU缓存。事实上，缓存将字节块从RAM加载到CPU寄存器。而后相邻元素就能高效地被加载了（顺序位置，或引用位置）。

3. 数据元素连续地存储在内存中，因此NumPy能够利用现代CPU的矢量化指令，像英特尔的SSE和AVX，AMD的XOP等。例如，为了做为CPU指令实现的矢量化算术计算，能够加载在128,256或512位寄存器中的多个连续的浮点数。

4. NumPy能够经过Intel Math Kernel Library (MKL)与高度优化的线性代数库相连，好比BLAS和LAPACK。NumPy中一些特定的矩阵计算也多是多线程，充分利用了现代多核处理器的优点。

总之，将数据存储在一个连续的内存块中，根据内存访问模式，CPU缓存和矢量化指令，能够确保以最佳方式使用现代CPU的体系结构。

就地操做和隐式拷贝操做之间的区别

让咱们解释一下技巧3。相似于a *= 2这样的表达式对应一个就地操做，即数组的全部元素值被乘以2。相比之下，a = a*2意味着建立了一个包含a*2结果值的新数组，变量a此时指向这个新数组。旧数组变为了无引用的，将被垃圾回收器删除。第一种状况中没有发生内存分配，相反，第二种状况中发生了内存分配。

更通常的状况，相似于a[i:j]这样的表达式是数组某些部分的视图：它们指向包含数据的内存缓冲区。利用就地操做改变它们，会改变原始数据。所以，a[:] = a * 2的结果是一个就地操做，和a = a * 2不同。

知道NumPy的这种细节能够帮助你解决一些错误（例如数组由于在一个视图上的一个操做，被无心中修改），并能经过减小没必要要的副本数量，优化代码的速度和内存消耗。

为何有些数组不进行拷贝操做，就不能被重塑？

一个转置的二维矩阵不依靠拷贝就没法进行铺平。一个二维矩阵包含的元素经过两个数字（行和列）进行索引，但它在内部是做为一个一维连续内存块存储的，可以使用一个数字访问。

有多个在一维内存块中存储矩阵元素的方法：咱们能够先放第一行的元素，而后第二行，以此类推，或者先放第一列的元素，而后第二列，以此类推。第一种方法叫作行优先排序，然后一种方法称为列优先排序。这两种方法之间的选择只是一个内部约定问题：NumPy使用行优先排序，相似于C，而不一样于FORTRAN。

更通常的状况，NumPy使用步长的概念进行多维索引和元素的底层序列（一维）内存位置之间的转换。array[i1, i2]和内部数据的相关字节地址之间的具体映射关系为：

1	offset = array.strides[0] * i1 + array.strides[1] * i2

重塑一个数组时，NumPy会尽量经过修改步长属性来避免拷贝。例如，当转置一个矩阵时，步长的顺序被翻转，但底层数据仍然是相同的。然而，仅简单地依靠修改步长没法完成铺平一个转置数组的操做（尝试下！），因此须要一个副本。

Recipe 4.6（NumPy中使用步长技巧）包含步长方面更普遍的讨论。同时，Recipe4.7（使用步长技巧实现一个高效的移动平均算法）展现了如何使用步伐加快特定数组计算。

from:http://blog.csdn.net/pipisorry/article/details/39087583

ref:Getting the Best Performance out of NumPy

IPython Cookbook

如何写出比 MATLAB 更快的矩阵运算程序？

Numpy使用MKL库提高计算性能