Python中优化NumPy包使用性能的教程
NumPy是Python中众多科学软件包的基础。它提供了一个特殊的数据类型ndarray，其在向量计算上做了优化。这个对象是科学数值计算中大多数算法的核心。
相比于原生的Python，利用NumPy数组可以获得显著的性能加速，尤其是当你的计算遵循单指令多数据流（SIMD）范式时。然而，利用NumPy也有可能有意无意地写出未优化的代码。
在这篇文章中，我们将看到一些技巧，这些技巧可以帮助你编写高效的NumPy代码。我们首先看一下如何避免不必要的数组拷贝，以节省时间和内存。因此，我们将需要深入NumPy的内部。
学习避免不必要的数据拷贝
NumPy数组计算可能涉及到内存块之间的内部拷贝。有时会有不必要的拷贝，此时应该避免。相应地这里有一些技巧，可以帮助你优化你的代码。
importnumpyasnp
查看数组的内存地址
1.查看静默数组拷贝的第一步是在内存中找到数组的地址。下边的函数就是做这个的：
defid(x): #Thisfunctionreturnsthememory #blockaddressofanarray. returnx.__array_interface__['data'][0]
2.有时你可能需要复制一个数组，例如你需要在操作一个数组时，内存中仍然保留其原始副本。
a=np.zeros(10);aid=id(a);aid 71211328 b=a.copy();id(b)==aid False
具有相同数据地址（比如id函数的返回值）的两个数组，共享底层数据缓冲区。然而，共享底层数据缓冲区的数组，只有当它们具有相同的偏移量（意味着它们的第一个元素相同）时，才具有相同的数据地址。共享数据缓冲区，但偏移量不同的两个数组，在内存地址上有细微的差别，正如下边的例子所展示的那样：
id(a),id(a[1:]) (71211328,71211336)
在这篇文章中，我们将确保函数用到的数组具有相同的偏移量。下边是一个判断两个数组是否共享相同数据的更可靠的方案：
defget_data_base(arr): """ForagivenNumpyarray,findsthebasearraythat"owns"theactualdata.""" base=arr whileisinstance(base.base,np.ndarray): base=base.base returnbase defarrays_share_data(x,y): returnget_data_base(x)isget_data_base(y) print(arrays_share_data(a,a.copy()),arrays_share_data(a,a[1:])) FalseTrue
感谢MichaelDroettboom指出这种更精确的方法，提出这个替代方案。
就地操作和隐式拷贝操作
3.数组计算包括就地操作（下面第一个例子：数组修改）或隐式拷贝操作（第二个例子：创建一个新的数组）。
a*=2;id(a)==aid True c=a*2;id(c)==aid False
一定要选择真正需要的操作类型。隐式拷贝操作很明显很慢，如下所示：
%%timeita=np.zeros(10000000) a*=2 10loops,bestof3:19.2msperloop %%timeita=np.zeros(10000000) b=a*2 10loops,bestof3:42.6msperloop
4.重塑一个数组可能涉及到拷贝操作，也可能涉及不到。原因将在下面解释。例如，重塑一个二维矩阵不涉及拷贝操作，除非它被转置（或更一般的非连续操作）：
a=np.zeros((10,10));aid=id(a);aid 53423728
重塑一个数组，同时保留其顺序，并不触发拷贝操作。
b=a.reshape((1,-1));id(b)==aid True
转置一个数组会改变其顺序，所以这种重塑会触发拷贝操作。
c=a.T.reshape((1,-1));id(c)==aid False
因此，后边的指令比前边的指令明显要慢。
5.数组的flatten和revel方法将数组变为一个一维向量（铺平数组）。flatten方法总是返回一个拷贝后的副本，而revel方法只有当有必要时才返回一个拷贝后的副本（所以该方法要快得多，尤其是在大数组上进行操作时）。
d=a.flatten();id(d)==aid False e=a.ravel();id(e)==aid True %timeita.flatten() 1000000loops,bestof3:881nsperloop %timeita.ravel() 1000000loops,bestof3:294nsperloop
广播规则
6.广播规则允许你在形状不同但却兼容的数组上进行计算。换句话说，你并不总是需要重塑或铺平数组，使它们的形状匹配。下面的例子说明了两个向量之间进行矢量积的两个方法：第一个方法涉及到数组的变形操作，第二个方法涉及到广播规则。显然第二个方法是要快得多。
n=1000 a=np.arange(n) ac=a[:,np.newaxis] ar=a[np.newaxis,:] %timeitnp.tile(ac,(1,n))*np.tile(ar,(n,1)) 100loops,bestof3:10msperloop %timeitar*ac 100loops,bestof3:2.36msperloop
在NumPy数组上进行高效的选择
NumPy提供了多种数组分片的方式。数组视图涉及到一个数组的原始数据缓冲区，但具有不同的偏移量，形状和步长。NumPy只允许等步长选择（即线性分隔索引）。NumPy还提供沿一个轴进行任意选择的特定功能。最后，花式索引（fancyindexing）是最一般的选择方法，但正如我们将要在文章中看到的那样，它同时也是最慢的。如果可能，我们应该选择更快的替代方法。
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)
4.比较一下两个方法的执行效率。
%timeita[::10] 1000000loops,bestof3:804nsperloop %timeita[np.arange(0,n,10)] 100loops,bestof3:14.1msperloop
花式索引慢好几个数量级，因为它要复制一个大数组。
替代花式索引：索引列表
5.当需要沿一个维度进行非等步长选择时，数组视图就无能为力了。然而，替代花式索引的方法在这种情况下依然存在。给定一个索引列表，NumPy的函数可以沿一个轴执行选择操作。
i=np.arange(0,n,10) b1=a[i] b2=np.take(a,i,axis=0) np.array_equal(b1,b2) True
第二个方法更快一点：
%timeita[i] 100loops,bestof3:13msperloop %timeitnp.take(a,i,axis=0) 100loops,bestof3:4.87msperloop
替代花式索引：布尔掩码
6.当沿一个轴进行选择的索引是通过一个布尔掩码向量指定时，compress函数可以作为花式索引的替代方案。
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 %timeita[i] 10loops,bestof3:59.8msperloop %timeitnp.compress(i,a,axis=0) 10loops,bestof3:24.1msperloop
第二个方法同样比花式索引快得多。
花式索引是进行数组任意选择的最一般方法。然而，往往会存在更有效、更快的方法，应尽可能首选那些方法。
当进行等步长选择时应该使用数组视图，但需要注意这样一个事实：视图涉及到原始数据缓冲区。
它是如何工作的？
在本节中，我们将看到使用NumPy时底层会发生什么，从而让我们理解该文章中的优化技巧。
为什么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位寄存器中的多个连续的浮点数。
此外，说一下这样一个事实：NumPy可以通过IntelMathKernelLibrary(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的这种细节可以帮助你解决一些错误（例如数组因为在一个视图上的一个操作，被无意中修改），并能通过减少不必要的副本数量，优化代码的速度和内存消耗。
为什么有些数组不进行拷贝操作，就不能被重塑？
我们在这里解释下技巧4，一个转置的二维矩阵不依靠拷贝就无法进行铺平。一个二维矩阵包含的元素通过两个数字（行和列）进行索引，但它在内部是作为一个一维连续内存块存储的，可使用一个数字访问。有多个在一维内存块中存储矩阵元素的方法：我们可以先放第一行的元素，然后第二行，以此类推，或者先放第一列的元素，然后第二列，以此类推。第一种方法叫做行优先排序，而后一种方法称为列优先排序。这两种方法之间的选择只是一个内部约定问题：NumPy使用行优先排序，类似于C，而不同于FORTRAN。
更一般的情况，NumPy使用步长的概念进行多维索引和元素的底层序列（一维）内存位置之间的转换。array[i1,i2]和内部数据的相关字节地址之间的具体映射关系为：
offset=array.strides[0]*i1+array.strides[1]*i2
重塑一个数组时，NumPy会尽可能通过修改步长属性来避免拷贝。例如，当转置一个矩阵时，步长的顺序被翻转，但底层数据仍然是相同的。然而，仅简单地依靠修改步长无法完成铺平一个转置数组的操作（尝试下！），所以需要一个副本。
Recipe4.6（NumPy中使用步长技巧）包含步长方面更广泛的讨论。同时，Recipe4.7（使用步长技巧实现一个高效的移动平均算法）展示了如何使用步伐加快特定数组计算。
内部数组排列还可以解释一些NumPy相似操作之间的意想不到的性能差异。作为一个小练习，你能解释一下下边这个例子吗？
a=np.random.rand(5000,5000) %timeita[0,:].sum() %timeita[:,0].sum() 100000loops,bestof3:9.57μsperloop 10000loops,bestof3:68.3μsperloop
NumPy的广播规则是什么？
广播规则描述了具有不同维度和/或形状的数组仍可以用于计算。一般的规则是：当两个维度相等，或其中一个为1时，它们是兼容的。NumPy使用这个规则，从后边的维数开始，向前推导，来比较两个元素级数组的形状。最小的维度在内部被自动延伸，从而匹配其他维度，但此操作并不涉及任何内存复制。