图解NumPy:常用函数的内在机制(上)

支持大量多维数组和矩阵运算的 NumPy 软件库是许多机器学习开发者和研究者的必备工具，本文将通过直观易懂的图示解析常用的 NumPy 功能和函数，帮助你理解 NumPy 操作数组的内在机制。

NumPy 是一个基础软件库，很多常用的 Python 数据处理软件库都使用了它或受到了它的启发，包括 pandas、PyTorch、TensorFlow、Keras 等。理解 NumPy 的工作机制能够帮助你提升在这些软件库方面的技能。而且在 GPU 上使用 NumPy 时，无需修改或仅需少量修改代码。

NumPy 的核心概念是 n 维数组。n 维数组的美丽之处是大多数运算看起来都一样，不管数组有多少维。但一维和二维有点特殊。本文分为三部分：

1. 向量：一维数组

2. 矩阵：二维数组

3. 三维及更高维

本文参考了 Jay Alammar 的文章《A Visual Intro to NumPy》并将其作为起点，然后进行了扩充，并做了一些细微修改。

NumPy 数组和 Python 列表

乍一看，NumPy 数组与 Python 列表类似。它们都可作为容器，能够快速获取和设置元素，但插入和移除元素会稍慢一些。

NumPy 数组完胜列表的最简单例子是算术运算：

除此之外，NumPy 数组的优势和特点还包括：

更紧凑，尤其是当维度大于一维时；

当运算可以向量化时，速度比列表更快；

当在后面附加元素时，速度比列表慢；

通常是同质的：当元素都是一种类型时速度很快。

这里 O(N) 的意思是完成该运算所需的时间和数组的大小成正比，而 O*(1)（即所谓的「均摊 O(1)」）的意思是完成运算的时间通常与数组的大小无关。

向量：一维数组

向量初始化

为了创建 NumPy 数组，一种方法是转换 Python 列表。NumPy 数组类型可以直接从列表元素类型推导得到。

要确保向其输入的列表是同一种类型，否则你最终会得到 dtype=’object’，这会影响速度，最终只留下 NumPy 中含有的语法糖。

NumPy 数组不能像 Python 列表一样增长。数组的末端没有留下任何便于快速附加元素的空间。因此，常见的做法是要么先使用 Python 列表，准备好之后再将其转换为 NumPy 数组，要么是使用 np.zeros 或 np.empty 预先留下必要的空间：

通常我们有必要创建在形状和元素类型上与已有数组匹配的空数组。

事实上，所有用于创建填充了常量值的数组的函数都带有 _like 的形式：

NumPy 中有两个函数能用单调序列执行数组初始化：

如果你需要类似 [0., 1., 2.] 这样的浮点数数组，你可以修改 arange 输出的类型：arange(3).astype(float)，但还有一种更好的方法。arange 函数对类型很敏感：如果你以整型数作为参数输入，它会生成整型数；如果你输入浮点数（比如 arange(3.)），它会生成浮点数。

但 arange 并不非常擅长处理浮点数：

在我们眼里，这个 0.1 看起来像是一个有限的十进制数，但计算机不这么看。在二进制表示下，0.1 是一个无限分数，因此必须进行约分，也由此必然会产生误差。也因为这个原因，如果向 arange 函数输入带分数部分的 step，通常得不到什么好结果：你可能会遇到差一错误 (off-by-one error)。你可以使该区间的末端落在一个非整数的 step 数中（solution1），但这会降低代码的可读性和可维护性。这时候，linspace 就可以派上用场了。它不受舍入的影响，总能生成你要求的元素数值。不过，使用 linspace 时会遇到一个常见的陷阱：它统计的是数据点的数量，而不是区间，因此其最后一个参数 num 通常比你所想的数大 1。因此，上面最后一个例子中的数是 11，而不是 10。

在进行测试时，我们通常需要生成随机数组：

向量索引

一旦你的数组中有了数据，NumPy 就能以非常巧妙的方式轻松地提供它们：

除了「花式索引（fancy indexing）」外，上面给出的所有索引方法都被称为「view」：它们并不存储数据，也不会在数据被索引后发生改变时反映原数组的变化情况。

所有包含花式索引的方法都是可变的：它们允许通过分配来修改原始数组的内容，如上所示。这一功能可通过将数组切分成不同部分来避免总是复制数组的习惯。

Python 列表与 NumPy 数组的对比

为了获取 NumPy 数组中的数据，另一种超级有用的方法是布尔索引（boolean indexing），它支持使用各类逻辑运算符：

any 和 all 的作用与在 Python 中类似，但不会短路。

不过要注意，这里不支持 Python 的「三元比较」，比如 3<=a<=5。

如上所示，布尔索引也是可写的。其两个常用功能都有各自的专用函数：过度重载的 np.where 函数和 np.clip 函数。它们的含义如下：

向量运算

NumPy 在速度上很出彩的一大应用领域是算术运算。向量运算符会被转换到 C++ 层面上执行，从而避免缓慢的 Python 循环的成本。NumPy 支持像操作普通的数那样操作整个数组。

与 Python 句法一样，a//b 表示 a 除 b（除法的商），x**n 表示 x?。

正如加减浮点数时整型数会被转换成浮点数一样，标量也会被转换成数组，这个过程在 NumPy 中被称为广播（broadcast）。

大多数数学函数都有用于处理向量的 NumPy 对应函数：

标量积有自己的运算符：

执行三角函数时也无需循环：

我们可以在整体上对数组进行舍入：

floor 为舍、ceil 为入，around 则是舍入到最近的整数（其中 .5 会被舍掉）

NumPy 也能执行基础的统计运算：

NumPy 的排序函数没有 Python 的排序函数那么强大：

Python 列表与 NumPy 数组的排序函数对比

在一维情况下，如果缺少 reversed 关键字，那么只需简单地对结果再执行反向，最终效果还是一样。二维的情况则会更困难一些（人们正在请求这一功能）。

搜索向量中的元素

与 Python 列表相反，NumPy 数组没有索引方法。人们很久之前就在请求这个功能，但一直还没实现。

Python 列表与 NumPy 数组的对比，index() 中的方括号表示可以省略 j 或同时省略 i 和 j。

一种查找元素的方法是 np.where(a==x)[0][0]，但这个方法既不优雅，速度也不快，因为它需要检查数组中的所有元素，即便所要找的目标就在数组起始位置也是如此。

另一种更快的方式是使用 Numba 来加速 next((i[0] for i, v in np.ndenumerate(a) if v==x), -1)。

一旦数组的排序完成，搜索就容易多了：v = np.searchsorted(a, x); return v if a[v]==x else -1 的速度很快，时间复杂度为 O(log N)，但它需要 O(N log N) 时间先排好序。

事实上，用 C 来实现它进而加速搜索并不是问题。问题是浮点比较。这对任何数据来说都不是一种简单直接可用的任务。

比较浮点数

函数 np.allclose(a, b) 能在一定公差下比较浮点数数组。

函数 np.allclose(a, b) 的工作过程示例。并没有万能方法！

np.allclose 假设所有被比较的数都在典型的 1 的范围内。举个例子，如果要在纳秒级的速度内完成计算，则需要用默认的 atol 参数值除以 1e9：np.allclose(1e-9, 2e-9, atol=1e-17) == False.

math.isclose 则不会对要比较的数进行任何假设，而是依赖用户给出合理的 abs_tol 值（对于典型的 1 的范围内的值，取默认的 np.allclose atol 值 1e-8 就足够好了）：math.isclose(0.1+0.2–0.3, abs_tol=1e-8)==True.

除此之外，np.allclose 在绝对值和相对公差的公式方面还有一些小问题，举个例子，对于给定的 a 和 b，存在 allclose(a, b) != allclose(b, a)。这些问题已在（标量）函数 math.isclose 中得到了解决，我们将在后面介绍它。对于这方面的更多内容，请参阅 GitHub 上的浮点数指南和对应的 NumPy 问题（https://floating-point-gui.de/errors/comparison/）。

矩阵：二维数组

NumPy 曾有一个专门的 matrix 类，但现在已经弃用了，所以本文会交替使用「矩阵」和「二维数组」这两个术语。

矩阵的初始化句法与向量类似：

这里必须使用双括号，因为第二个位置参数是 dtype（可选，也接受整数）。

随机矩阵生成的句法也与向量的类似：

二维索引的句法比嵌套列表更方便：

view 符号的意思是当切分一个数组时实际上没有执行复制。当该数组被修改时，这些改变也会反映到切分得到的结果上。

axis 参数

在很多运算中（比如 sum），你需要告诉 NumPy 是在列上还是行上执行运算。为了获取适用于任意维度的通用符号，NumPy 引入了 axis 的概念：事实上，axis 参数的值是相关问题中索引的数量：第一个索引为 axis=0，第二个索引为 axis=1，以此类推。因此在二维情况下，axis=0 是按列计算，axis=1 是按行计算。

矩阵算术运算

除了逐元素执行的常规运算符（比如 +、-、、/、//、*），这里还有一个计算矩阵乘积的 @ 运算符：

我们已在第一部分介绍过标量到数组的广播，在其基础上进行泛化后，NumPy 支持向量和矩阵的混合运算，甚至两个向量之间的运算：

二维数组中的广播

行向量和列向量

正如上面的例子所示，在二维情况下，行向量和列向量的处理方式有所不同。这与具备某类一维数组的 NumPy 实践不同（比如二维数组 a— 的第 j 列 a[:,j] 是一个一维数组）。默认情况下，一维数组会被视为二维运算中的行向量，因此当用一个矩阵乘以一个行向量时，你可以使用形状 (n,) 或 (1, n)——结果是一样的。如果你需要一个列向量，则有多种方法可以基于一维数组得到它，但出人意料的是「转置」不是其中之一。

基于一维数组得到二维数组的运算有两种：使用 reshape 调整形状和使用 newaxis 进行索引：

其中 -1 这个参数是告诉 reshape 自动计算其中一个维度大小，方括号中的 None 是用作 np.newaxis 的快捷方式，这会在指定位置添加一个空 axis。

因此，NumPy 共有三类向量：一维向量、二维行向量和二维列向量。下图展示了这三种向量之间的转换方式：

一维向量、二维行向量和二维列向量之间的转换方式。根据广播的原则，一维数组可被隐含地视为二维行向量，因此通常没必要在这两者之间执行转换——因此相应的区域被阴影化处理。