了解NumPy的内部原理，避免不必要的数组复制
来源于： IPython Cookbook, Second Edition, by Cyrille Rossant
▶ Code on GitHub with a MIT license
▶ Go to Chapter 4 : Profiling and Optimization
▶ Get the Jupyter notebook
前言
近期再搞大数据相关的事情， 发现大数据，最大的问题是大 数据，数据量大然后就会出现各种问题， 之前提到的空间和效率问题， 这也是属于NP内部的空间和效率的一种问题，之前没有太多的了解，一头扎进去，发现不明白其中的内部原理才是最大的问题， 不明白数据结构和数据存储处理过程，走了很多弯路，这里记录下。
与本机Python代码相比，可以使用NumPy实现显著的性能加速，特别是当计算遵循单指令多数据(SIMD)范式时。然而，也有可能无意中使用NumPy编写未优化的代码。
在接下来的几个case中，将看到一些技巧，可以帮助编写优化的NumPy代码。在本case中，将了解如何避免不必要的数组复制以节省内存。在这方面，需要深入了解NumPy的内部结构。
准备阶段
首先，需要一种方法来检查两个数组在内存中是否共享相同的底层数据缓冲区。定义一个函数aid()，它返回底层数据缓冲区的内存位置:
import numpy as np
def aid(x):
# This function returns the memory
# block address of an array.
return x.__array_interface__['data'][0]
具有相同数据位置(由aid()返回)的两个数组共享相同的底层数据缓冲区，只有当数组具有相同的偏移量(意味着它们具有相同的第一个元素)时。具有不同偏移量的两个共享数组的内存位置略有不同，如下例所示:
a = np.zeros(3)
aid(a), aid(a[1:])
(21535472, 21535480)
在接下来的几个case中，将确保对具有相同偏移量的数组使用此方法。下面是一个更通用、更可靠的解决方案，用于确定两个数组是否共享相同的数据:
def get_data_base(arr):
"""For a given NumPy array, find the base array
that owns the actual data."""
base = arr
while isinstance(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()))
False
print(arrays_share_data(a, a[:1]))
True
怎么验证
NumPy数组的计算可能涉及内存块之间的内部拷贝。这些拷贝并不总是必要的，在这种情况下应该避免，正如将在以下提示中看到的那样。
1. 有时可能需要复制一个数组;例如，如果需要操作一个数组，同时在内存中保留一个原始副本:
import numpy as np
a = np.zeros(10)
ax = aid(a)
ax
32250112
b = a.copy()
aid(b) == ax
False
- 数组计算可以包括原地操作(以下代码中的第一个示例:数组被修改)或隐式复制操作(第二个示例:创建一个新数组):
a *= 2
aid(a) == ax
True
c = a * 2
aid(c) == ax
False
隐式复制操作比较慢，如下所示:
%%timeit a = np.zeros(10000000)
a *= 2
4.85 ms ± 24 µs per loop (mean ± std. dev. of 7 runs,
100 loops each)
%%timeit a = np.zeros(10000000)
b = a * 2
7.7 ms ± 105 µs per loop (mean ± std. dev. of 7 runs,
100 loops each)
3.重新塑造数组可能涉及，也可能不涉及副本。原因将在它的工作原理中解释这部分。例如，重塑一个2D矩阵不涉及复制，除非它是转置(或更一般地说，不连续):
#首先重塑
a = np.zeros((100, 100))
ax = aid(a)
b = a.reshape((1, -1))
aid(b) == ax
True
#转秩
c = a.T.reshape((1, -1))
aid(c) == ax
False
因此，后者指令明显比前者慢:
%timeit b = a.reshape((1, -1))
330 ns ± 0.517 ns per loop (mean ± std. dev. of 7 runs
1000000 loops each)
%timeit a.T.reshape((1, -1))
5 µs ± 5.68 ns per loop (mean ± std. dev. of 7 runs,
100000 loops each)
4. 数组的flatten()和ravel()方法都将其重塑为一维矢量(一个扁平数组)。然而，flatten()方法总是返回一个副本，而ravel()方法只在必要时返回一个副本(因此它更快，特别是对于大型数组)。
d = a.flatten()
aid(d) == ax
False
e = a.ravel()
aid(e) == ax
True
%timeit a.flatten()
2.3 µs ± 18.1 ns per loop (mean ± std. dev. of 7 runs,
100000 loops each)
%timeit a.ravel()
199 ns ± 5.02 ns per loop (mean ± std. dev. of 7 runs,
10000000 loops each)
5. 广播规则允许在具有不同但兼容形状的数组上进行计算。换句话说，并不总是需要重塑或平铺数组以使它们的形状匹配。下面的例子说明了在两个向量之间做外积的两种方法:第一种方法涉及到数组平铺，第二种(更快)涉及到广播:
n = 1000
a = np.arange(n)
ac = a[:, np.newaxis] # column vector
ar = a[np.newaxis, :] # row vector
%timeit np.tile(ac, (1, n)) * np.tile(ar, (n, 1))
5.7 ms ± 42.6 µs per loop (mean ± std. dev. of 7 runs,
100 loops each)
%timeit ar * ac
784 µs ± 2.39 µs per loop (mean ± std. dev. of 7 runs,
1000 loops each)
为什么会这样
在本节中，将了解在使用NumPy时发生了什么，以及这些知识如何让我们理解这个case中给出的技巧。
为什么NumPy数组是高效的?
NumPy数组基本上是由元数据(特别是维数、形状和数据类型)和实际数据描述的。数据存储在一个均匀且连续的内存块中，位于系统内存(随机存取内存或RAM)中的特定地址。这块内存称为数据缓冲区。这是数组和纯Python结构(如列表)之间的主要区别，列表中的项分散在系统内存中。这方面是NumPy数组如此高效的关键特性。
为什么这一点如此重要?以下是主要原因:
- 数组上的计算可以用低级语言(比如C语言)非常有效地编写(NumPy的很大一部分实际上是用C语言编写的)。例如，知道内存块的地址和数据类型，就可以简单地遍历所有项。在Python中使用列表执行此操作将会有很大的开销。
- 由于CPU缓存，内存访问模式中的空间局部性可以显著提高性能。事实上，缓存将字节块从RAM加载到CPU寄存器。然后非常有效地加载相邻项(顺序局部性或引用局部性)。
- 最后，项目连续存储在内存中的事实允许NumPy利用现代cpu的向量化指令，如英特尔的SSE和AVX, AMD的XOP等。例如，多个连续浮点数可以装入128位、256位或512位寄存器，用于作为CPU指令实现的向量化算术计算。
此外，NumPy可以通过ATLAS或英特尔数学内核库(MKL)链接到高度优化的线性代数库，如BLAS和LAPACK。一些特定的矩阵计算也可以是多线程的，利用现代多核处理器的强大功能。
总之，将数据存储在一个连续的内存块中可以确保现代CPU的架构在内存访问模式、CPU缓存和向量化指令方面得到最佳使用。
in-place就地复制操作和implicit-copy隐式复制操作之间的区别是什么?
让我们在第2步中解释这个例子。像a *= 2这样的表达式对应于就地操作，其中数组的所有值都乘以2。相比之下，a = a*2意味着创建了一个包含a*2值的新数组，变量a现在指向这个新数组。旧数组不再被引用，将被垃圾回收器删除。与第二种情况相反，第一种情况下没有发生内存分配。
更一般地说，像a[i:j]这样的表达式是对数组部分的视图;它们指向包含数据的内存缓冲区。用就地操作修改它们会改变原始数组。
了解NumPy的这种微妙之处可以帮助您修复一些错误(由于视图上的操作而隐式地无意地修改了数组)，并通过减少不必要的副本数量来优化代码的速度和内存消耗。
为什么一些数组不能在没有副本的情况下被重新塑造?
在这里的第3步中解释了一个例子，其中一个转置的2D矩阵不能在没有副本的情况下被平坦化。2D矩阵包含以两个数字(行和列)为索引的项，但它在内部存储为一个1D连续内存块，仅用一个数字访问。在一维内存块中存储矩阵项的方法不止一种:可以将第一行的元素放在前面，然后是第二行，等等，或者将第一列的元素放在前面，然后是第二列，等等。第一种方法称为行-主序，而后者称为列-主序。在这两种方法之间进行选择只是内部约定的问题:NumPy像C一样使用行主序，但与FORTRAN不同。
通常，NumPy使用跨行的概念在多维索引和底层(1D)元素序列的内存位置之间进行转换。数组[i1, i2]与内部数据的相关字节地址之间的具体映射由:
offset = array.strides[0] * i1 + array.strides[1] * i2
在重新塑造数组时，NumPy通过修改strides属性尽可能避免复制。例如，当转置一个矩阵时，步长顺序颠倒，但底层数据保持相同。然而，仅仅通过修改步长是无法实现转置数组的扁平化的，因此需要一个副本。
内部数组布局还可以解释非常相似的NumPy操作之间的一些意想不到的性能差异。作为一个小练习，能解释一下下面的基准测试吗?
a = np.random.rand(5000, 5000)
%timeit a[0, :].sum()
2.91 µs ± 20 ns per loop (mean ± std. dev. of 7 runs,
100000 loops each)
%timeit a[:, 0].sum()
33.7 µs ± 22.7 ns per loop (mean ± std. dev. of 7 runs
10000 loops each)
简单解释下， 已经存在的array a ， a[0,:].sum()相当于遍历一行求和， a[:, 0].sum()相当于遍历一列求和，行数据在内存中挨着存储，列数据在内存中是间隔存储的。
什么是NumPy广播规则?
广播规则描述了如何使用不同尺寸和/或形状的数组进行计算。一般的规则是，当两个维度相等时，或者当其中一个维度为1时，它们是兼容的。NumPy使用这个规则来比较两个数组元素的形状，从后面的维度开始，一直向前。最小的维度在内部拉伸以匹配其他维度，但此操作不涉及任何内存复制。
以下是一些参考资料:
广播规则和示例：
http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html
NumPy中的数组接口:
http://docs.scipy.org/doc/numpy/reference/arrays.interface.html
参考地点
https://en.wikipedia.org/wiki/Locality_of_reference
SciPy讲座笔记中的NumPy内部知识
http://scipy-lectures.github.io/advanced/advanced_numpy
Nicolas Rougier的100个NumPy练习
http://www.loria.fr/~rougier/teaching/numpy.100/index.html