NumPy数值计算#

计算机毫无用处,它们只能给出答案。

—Pablo Picasso

虽然Python解释器本身已经提供了丰富的数据结构,但NumPy和其他库以有价值的方式进一步扩展了这些数据结构。本章聚焦于NumPy,它提供了多维数组对象(Multidimensional Array)来存储同质或异质数据数组,并支持代码的向量化。

本章涵盖以下数据结构:

对象类型含义用途
ndarray(常规)n维数组对象大型数值数据数组
ndarray(记录)二维数组对象按列组织的表格数据

本章组织如下:

"数据数组" on page 86

本节讨论使用纯Python代码处理数据数组。

"常规NumPy数组" on page 90

这是关于常规NumPy ndarray类的核心部分,是所有涉及数值数据的数据密集型Python用例的主力工具。

"结构化NumPy数组" on page 105

本小节简要介绍用于处理带列名的表格数据的结构化(或记录)ndarray对象。

"代码向量化" on page 106

本节讨论代码向量化及其优势;还讨论了某些场景下内存布局的重要性。

数据数组#

前一章展示了Python提供了一些相当有用和灵活的通用数据结构。特别是,列表对象可以被视为具有许多便利特性和应用领域的真正主力。使用这种灵活(可变)的数据结构是有代价的,表现为相对较高的内存使用、较慢的性能,或两者兼有。然而,科学和金融应用通常需要在特殊数据结构上进行高性能操作。在这方面最重要的数据结构之一就是数组(Array)。数组通常将相同数据类型的其他(基本)对象按行和列组织起来。

暂时假设只有数字是相关的,尽管这个概念也可以推广到其他类型的数据。在最简单的情况下,一维数组在数学上表示一个向量,通常由浮点数对象在内部表示的真实数组成。它仅由单行或单列元素组成。在更常见的情况下,数组表示一个 i × j 的元素矩阵。这个概念可以推广到三维中的 i × j × k 立方体元素,以及一般的 n 维数组,其形状为 $i \times j \times k \times l \times$

线性代数和向量空间理论等数学学科表明,这种数学结构在许多科学学科和领域中具有高度重要性。因此,拥有一个专门设计用于方便高效处理数组的数据结构类是富有成效的。这正是Python库NumPy及其强大的 ndarray 类的用武之地。在下一节介绍这个类之前,本节将说明处理数组的两种替代方案。

使用Python列表的数组#

可以使用前一章介绍的内置数据结构来构建数组。list 对象特别适合完成此任务。一个简单的列表已经可以被视为一维数组:

In [1]: v = [0.5, 0.75, 1.0, 1.5, 2.0]

包含数字的列表对象。

由于列表对象可以包含任意其他对象,它们也可以包含其他列表对象。通过这种方式,可以通过嵌套列表对象轻松构造二维和更高维的数组:

In [5]: v1 = [0.5, 1.5]
    v2 = [1, 2]
    m = [v1, v2]
    c = [m, m] ①
    c
Out[5]: [[0.5, 1.5], [1, 2]], [[0.5, 1.5], [1, 2]])
In [6]: c[1][1][0]
Out[6]: 1
In [2]: m = [v, v, v] ①
    m ②
Out[2]: [[0.5, 0.75, 1.0, 1.5, 2.0],
    [0.5, 0.75, 1.0, 1.5, 2.0],
    [0.5, 0.75, 1.0, 1.5, 2.0]]

包含列表对象的列表对象……

……从而得到一个数字矩阵。

还可以通过简单索引轻松选择行,或通过双重索引选择单个元素(但整个列并不容易选择):

In [3]: m[1]
Out[3]: [0.5, 0.75, 1.0, 1.5, 2.0]
In [4]: m[1][0]
Out[4]: 0.5

嵌套可以进一步推广到更通用的结构:

数字立方体。

注意,以刚才的方式组合对象通常使用的是指向原始对象的引用指针。这在实践中意味着什么?看一下以下操作:

In [7]: v = [0.5, 0.75, 1.0, 1.5, 2.0]
    m = [v, v, v]
    m
Out[7]: [[0.5, 0.75, 1.0, 1.5, 2.0],
    [0.5, 0.75, 1.0, 1.5, 2.0],
    [0.5, 0.75, 1.0, 1.5, 2.0]]

现在更改 v 对象第一个元素的值,看看 m 对象会发生什么:

In [8]: v[0] = 'Python'
m
Out[8]: [['Python', 0.75, 1.0, 1.5, 2.0],
['Python', 0.75, 1.0, 1.5, 2.0],
['Python', 0.75, 1.0, 1.5, 2.0]]
In [9]: from copy import deepcopy
    v = [0.5, 0.75, 1.0, 1.5, 2.0]
    m = 3 * [deepcopy(v), ]
    m
Out[9]: [[0.5, 0.75, 1.0, 1.5, 2.0],
    [0.5, 0.75, 1.0, 1.5, 2.0],
    [0.5, 0.75, 1.0, 1.5, 2.0]]
In [10]: v[0] = 'Python'
    m
Out[10]: [[0.5, 0.75, 1.0, 1.5, 2.0],
    [0.5, 0.75, 1.0, 1.5, 2.0],
    [0.5, 0.75, 1.0, 1.5, 2.0]]

可以通过使用 copy 模块的 deepcopy() 函数来避免这种情况:

使用物理副本而非引用指针。

因此,原始对象的更改……

……不再有任何影响。

Python的array类#

Python中有一个专门的 array 模块。根据文档:

该模块定义了一种对象类型,可以紧凑地表示基本值数组:字符、整数、浮点数。数组是序列类型,行为与列表非常相似,只是存储的对象类型受到约束。类型在对象创建时通过类型码指定,类型码是一个单个字符。

考虑以下代码,它从一个列表对象实例化一个数组对象:

In [11]: v = [0.5, 0.75, 1.0, 1.5, 2.0]
In [12]: import array
In [13]: a = array.array('f', v) ①
a
Out[13]: array('f', [0.5, 0.75, 1.0, 1.5, 2.0])
In [14]: a.append(0.5) ②
a
Out[14]: array('f', [0.5, 0.75, 1.0, 1.5, 2.0, 0.5])
In [15]: a.extend([5.0, 6.75]) ②
a
Out[15]: array('f', [0.5, 0.75, 1.0, 1.5, 2.0, 0.5, 5.0, 6.75])
In [16]: 2 * a ③
Out[16]: array('f', [0.5, 0.75, 1.0, 1.5, 2.0, 0.5, 5.0, 6.75, 0.5, 0.75, 1.0, 1.5, 2.0, 0.5, 5.0, 6.75])

使用 float 作为类型码实例化数组对象。

主要方法与列表对象的方法类似。

虽然"标量乘法"原则上可行,但结果不是数学上预期的结果;相反,元素被重复了。

尝试追加与指定数据类型不同类的对象会引发 TypeError:

In [17]: a.append('string') ①
TypeErrorTraceback (most recent call last)
<ipython-input-17-14cd6281866b> in <module>()
----> 1 a.append('string') ①
TypeError: must be real number, not str
In [18]: a.tolist() ②
Out[18]: [0.5, 0.75, 1.0, 1.5, 2.0, 0.5, 5.0, 6.75]

只能追加 float 对象;其他数据类型/类型码会引发错误。

然而,如果需要这样的灵活性,数组对象可以轻松转换回列表对象。

array 类的一个优点是它具有内置的存储和检索功能:

In [19]: f = open('array.apy', 'wb') ①
    a-.tofile(f) ②
    f.close() ③

In [20]: with open('array.apy', 'wb') as f: ④
    a-.tofile(f) ④

In [21]: !ls -n arr* ⑤
    -rw-r--r--@ 15032032 Nov 711:46 array.apy

打开一个文件用于写入二进制数据。

将数组数据写入文件。

关闭文件。

备选方案:使用 with 上下文执行相同操作。

显示写入磁盘的文件。

如前所述,从磁盘读取数据时,数组对象的数据类型也很重要:

In [22]: b = array.array('f') ①
In [23]: with open('array.apy', 'rb') as f: ②
    b.fromfile(f, 5) ③
In [24]: b ③
Out[24]: array('f', [0.5, 0.75, 1.0, 1.5, 2.0])
In [25]: b = array.array('d') ④
In [26]: with open('array.apy', 'rb') as f:
    b.fromfile(f, 2) ⑤
In [27]: b ⑥
Out[27]: array('d', [0.0004882813645963324, 0.12500002956949174])

实例化一个类型码为 float 的新数组对象。

以读取二进制模式打开文件……

……并将五个元素读入 b 对象。

实例化一个类型码为 double 的新数组对象。

从文件中读取两个元素。

类型码的差异导致"错误"的数字。

常规NumPy数组#

用列表对象组合数组结构是可行的,但并不是很方便,而且 list 类的设计初衷也不是为了实现这一目标。它的范围更广泛、更通用。array 类稍专门化一些,提供了一些处理数据数组的有用特性。然而,一个真正专门化的类对于处理数组类型的结构将非常有益。

基础知识#

numpy.ndarray 正是这样一个类,其设计目标便是方便、高效地处理 n 维数组——即以高性能的方式。通过示例最能说明该类实例的基本使用方法:

In [28]: import numpy as np ①
In [29]: a = np.array([0, 0.5, 1.0, 1.5, 2.0]) ②
a
Out[29]: array([0., 0.5, 1., 1.5, 2.])
In [30]: type(a) ②
Out[30]: numpy.ndarray
In [31]: a = np.array(['a', 'b', 'c']) ③
a
Out[31]: array(['a', 'b', 'c'], dtype='<U1')
In [32]: a = np.arange(2, 20, 2) ④
a
Out[32]: array([2, 4, 6, 8, 10, 12, 14, 16, 18])
In [33]: a = np.arange(8, dtype=np.float) ⑤
a
Out[33]: array([0., 1., 2., 3., 4., 5., 6., 7.])
In [34]: a[5:] ⑥
Out[34]: array([5., 6., 7.])
In [35]: a[:2] ⑥
Out[35]: array([0., 1.])

导入 numpy 包。

从一个包含浮点数的列表对象创建 ndarray 对象。

从一个包含字符串的列表对象创建 ndarray 对象。

np.arange() 类似于 range()……

……但额外接受 dtype 参数。

对于一维 ndarray 对象,索引操作与通常相同。

ndarray 类的一个主要特点是其众多的内置方法。例如:

In [36]: a.sum() ①
Out[36]: 28.0
In [37]: a.std() ②
Out[37]: 2.29128784747792
In [38]: a.cumsum() ③
Out[38]: array([0., 1., 3., 6., 10., 15., 21., 28.])

所有元素的和。

元素的标准差。

所有元素的累积和(从索引位置0开始)。

另一个主要特性是定义在 ndarray 对象上的(向量化)数学运算:

In [39]: l = [0., 0.5, 1.5, 3., 5.]
    2 * l
Out[39]: [0.0, 0.5, 1.5, 3.0, 5.0, 0.0, 0.5, 1.5, 3.0, 5.0]

In [40]: a
Out[40]: array([0., 1., 2., 3., 4., 5., 6., 7.])

In [41]: 2 * a
Out[41]: array([0., 2., 4., 6., 8., 10., 12., 14.])

In [42]: a ** 2
Out[42]: array([0., 1., 4., 9., 16., 25., 36., 49.])

In [43]: 2 ** a
Out[43]: array([1., 2., 4., 8., 16., 32., 64., 128.])

In [44]: a ** a
Out[44]: array([1.00000e+00, 1.00000e+00, 4.00000e+00, 2.70000e+01, 2.56000e+02, 3.12500e+03, 4.66560e+04, 8.23543e+05])

列表对象的标量乘法导致元素重复。

相比之下,ndarray 对象实现了真正的标量乘法。

这会逐元素计算平方值。

这将 ndarray 的元素解释为指数。

这计算每个元素的自乘幂。

通用函数(Universal Functions)是 NumPy 包的另一个重要特性。它们之所以被称作"通用",是因为通常既可以操作 ndarray 对象,也可以操作基本Python数据类型。然而,当将通用函数应用于(例如)Python float 对象时,需要注意性能比 math 模块中的相同功能要低:

In [50]: math.sqrt(a) 5
In [45]: np.exp(a)
Out[45]: array([1.00000000e+00, 2.71828183e+00, 7.38905610e+00, 2.00855369e+01, 5.45981500e+01, 1.48413159e+02, 4.03428793e+02, 1.09663316e+03])
In [46]: np.sqrt(a) ②
Out[46]: array([0., 1., 1.41421356, 1.73205081, 2.2.23606798, 2.44948974, 2.64575131])
In [47]: np.sqrt(2.5) ③
Out[47]: 1.5811388300841898
In [49]: math.sqrt(2.5) 4
Out[49]: 1.5811388300841898
TypeErrorTraceback (most recent call last)
<ipython-input-50-b39de4150838> in <module>()
----> 1 math.sqrt(a) ⑤

TypeError: only size-1 arrays can be converted to Python scalars

In [51]: %timeit np.sqrt(2.5) 722 ns ± 13.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

In [52]: %timeit math.sqrt(2.5) ⑦
91.8 ns ± 4.13 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)

逐元素计算指数值。

为每个元素计算平方根。

计算Python float 对象的平方根。

同样的计算,这次使用 math 模块。

math.sqrt() 函数不能直接应用于 ndarray 对象。

将通用函数 np.sqrt() 应用于Python float 对象……

……比使用 math.sqrt() 函数的相同操作慢得多。

多维#

过渡到多个维度是无缝的,到目前为止介绍的所有特性都可以推广到更一般的情况。特别是,索引系统在所有维度上保持一致:

In [53]: b = np.array([a, a * 2]) ①
b
Out[53]: array([[0., 1., 2., 3., 4., 5., 6., 7.], [0., 2., 4., 6., 8., 10., 12., 14.]])

In [54]: b[0] ②
Out[54]: array([0., 1., 2., 3., 4., 5., 6., 7.])

In [55]: b[0, 2] ③
Out[55]: 2.0

In [56]: b[:, 1] ④
Out[56]: array([1., 2.])

In [57]: b.sum() ⑤
Out[57]: 84.0

In [58]: b.sum(axis=0) ⑥
Out[58]: array([0., 3., 6., 9., 12., 15., 18., 21].])

In [59]: b.sum(axis=1) ⑦
Out[59]: array([28., 56.])

从一维 ndarray 对象构造二维 ndarray 对象。

选择第一行。

选择第一行中的第三个元素;索引在方括号内用逗号分隔。

选择第二列。

计算所有值的总和。

沿第一个轴计算总和,即按列求和。

沿第二个轴计算总和,即按行求和。

有多种方法可以初始化(实例化)ndarray 对象。一种是如前所述,通过 np.array。但这假设数组的所有元素都已经可用。相反,可能希望先实例化 ndarray 对象,然后再用代码执行过程中生成的结果来填充它们。为此,可以使用以下函数:

In [60]: c = np.zeros((2, 3), dtype='i', order='C') ①
    c
Out[60]: array([[0, 0, 0], [0, 0, 0]], dtype=int32)

In [61]: c = np.ones((2, 3, 4), dtype='i', order='C') ②
    c
Out[61]: array([[1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]],
    [[1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]], dtype=int32)

In [62]: d = np.zeros_like(c, dtype='f16', order='C') ③
    d
Out[62]: array([[0., 0., 0., 0.], [0., 0., 0., 0.], [0., 0., 0., 0.]],
    [[0., 0., 0., 0.], [0., 0., 0., 0.], [0., 0., 0., 0.]]], dtype=float128)

In [63]: d = np.ones_like(c, dtype='f16', order='C') ③
    d
Out[63]: array([[1., 1., 1., 1.], [1., 1., 1., 1.], [1., 1., 1., 1.]],
    [[1., 1., 1., 1.], [1., 1., 1., 1.], [1., 1., 1., 1.]]], dtype=float128)

In [64]: e = np.empty((2, 3, 2)) ④
    e
Out[64]: array([[0.00000000e+000, 0.00000000e+000], [0.00000000e+000, 0.00000000e+000], [0.00000000e+000, 0.00000000e+000]],
    [[0.0000000e+000, 0.0000000e+000], [0.0000000e+000, 7.49874326e+247], [1.28822975e-231, 4.33190018e-311]])]

In [65]: f = np.empty_like(c) ④
    f
Out[65]: array([[    0,    0,    0,
    [    0,    0,    0,
[ 0, 0, 0, 0]],
[[ 0, 0, 0, 0],
[ 0, 0, 740455269, 1936028450],
[ 0, 268435456, 1835316017, 2041]]], dtype=int32)
In [66]: np.eye(5) ⑤
Out[66]: array([[1., 0., 0., 0., 0.],
[0., 1., 0., 0., 0.],
[0., 0., 1., 0., 0.],
[0., 0., 0., 1., 0.],
[0., 0., 0., 0., 1.]])
In [67]: g = np.linspace(5, 15, 12) ⑥
g
Out[67]: array([ 5. , 5.90909091, 6.81818182, 7.72727273, 8.63636364,
9.54545455, 10.45454545, 11.36363636, 12.27272727, 13.18181818,
14.09090909, 15. ])

创建预填充为零的 ndarray 对象。

创建预填充为一的 ndarray 对象。

类似功能,但从另一个 ndarray 对象推断形状。

创建不预填充任何内容的 ndarray 对象(数字取决于内存中的位模式)。

创建一个对角线填充为一的方阵,作为 ndarray 对象。

创建一维 ndarray 对象,元素之间均匀间隔;参数为起始值、结束值和数量。

对于所有这些函数,可以提供以下参数:

shape#

一个 int、int 对象的序列,或对另一个 ndarray 的引用

dtype(可选)#

dtype——这些是 NumPy 特定的 ndarray 对象数据类型

order(可选)#

元素在内存中的存储顺序:C 表示类似 C 语言(即行优先)或 F 表示类似 Fortran(即列优先)

在这里,与基于 list 的方法相比,NumPy 如何使用 ndarray 类专门化数组的构建变得显而易见:

• ndarray 对象具有内置的维度(轴)。

• ndarray 对象是不可变的;其长度(大小)是固定的。

• 它只允许整个数组使用单一数据类型(np.dtype)。

相比之下,array 类仅共享允许单一数据类型(类型码、dtype)的特性。

order 参数的作用将在本章后面讨论。表4-1概述了选定的 np.dtype 对象(即 NumPy 允许的基本数据类型)。

表4-1. NumPy dtype 对象

dtype描述示例
?布尔值? (True 或 False)
i有符号整数i8 (64位)
u无符号整数u8 (64位)
f浮点数f8 (64位)
c复数浮点数c32 (256位)
m时间差m (64位)
M日期时间M (64位)
O对象O (指向对象的指针)
UUnicodeU24 (24个Unicode字符)
V原始数据(空类型)V12 (12字节数据块)

元信息#

每个 ndarray 对象都提供对一些有用属性的访问:

In [68]: g.size ①
Out[68]: 12

In [69]: g.itemsize ②
Out[69]: 8

In [70]: g.ndim ③
Out[70]: 1

In [71]: g.shape ④
Out[71]: (12,)
In [72]: g.dtype ⑤
Out[72]: dtype('float64')

In [73]: g.nbytes ⑥
Out[73]: 96

元素数量。

表示一个元素所使用的字节数。

维度数量。

ndarray 对象的形状。

元素的 dtype。

内存中使用的总字节数。

重塑与大小调整#

虽然 ndarray 对象默认是不可变的,但有多种选项可以重塑(Reshape)和调整大小(Resize)这样的对象。重塑通常只是提供对同一数据的另一种视图,而调整大小通常会创建一个新的(临时)对象。首先,一些重塑的例子:

In [74]: g = np.arange(15)

In [75]: g
Out[75]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])

In [76]: g.shape ①
Out[76]: (15,)

In [77]: np.shape(g) ①
Out[77]: (15,)

In [78]: g.reshape((3, 5)) ②
Out[78]: array([[0, 1, 2, 3, 4],
    [5, 6, 7, 8, 9],
    [10, 11, 12, 13, 14]])
In [79]: h = g.reshape((5, 3)) ③
h
Out[79]: array([[0, 1, 2],
    [3, 4, 5],
    [6, 7, 8],
    [9, 10, 11],
    [12, 13, 14]])
In [80]: h.T ④
Out[80]: array([[0, 3, 6, 9, 12],
    [1, 4, 7, 10, 13],
    [2, 5, 8, 11, 14]])
In [81]: h.transpose() ④
Out[81]: array([[0, 3, 6, 9, 12],
    [1, 4, 7, 10, 13],
    [2, 5, 8, 11, 14]])

原始 ndarray 对象的形状。

重塑为二维(内存视图)。

创建一个新对象。

新 ndarray 对象的转置。

在重塑操作期间,ndarray 对象中的元素总数不变。在调整大小操作期间,这个数量会改变——要么减少(缩小)要么增加(放大)。以下是一些调整大小的例子:

In [82]: g
Out[82]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])

In [83]: np.resize(g, (3, 1)) ①
Out[83]: array([[0],
    [1],
    [2]])

In [84]: np.resize(g, (1, 5)) ①
Out[84]: array([[0, 1, 2, 3, 4]])
In [85]: np.resize(g, (2, 5)) ①
Out[85]: array([[0, 1, 2, 3, 4],
    [5, 6, 7, 8, 9]])
In [86]: n = np.resize(g, (5, 4)) ②
n
Out[86]: array([[0, 1, 2, 3],
    [4, 5, 6, 7],
    [8, 9, 10, 11],
    [12, 13, 14, 0],
    [1, 2, 3, 4]])

二维,缩小。

二维,放大。#

堆叠(Stacking)是一种特殊操作,允许水平或垂直组合两个 ndarray 对象。但是,"连接"维度的大小必须相同:

In [87]: h
Out[87]: array([[0, 1, 2], [3, 4, 5],
[ 6, 7, 8],
[ 9, 10, 11],
[12, 13, 14]])

In [88]: np.hstack((h, 2 * h)) ①
Out[88]: array([[ 0, 1, 2, 0, 2, 4],
    [ 3, 4, 5, 6, 8, 10],
    [ 6, 7, 8, 12, 14, 16],
    [ 9, 10, 11, 18, 20, 22],
    [12, 13, 14, 24, 26, 28]])

In [89]: np.vstack((h, 0.5 * h)) ②
Out[89]: array([[ 0., 1., 2. ],
    [ 3., 4., 5. ],
    [ 6., 7., 8. ],
    [ 9., 10., 11. ],
    [12., 13., 14. ],
    [ 0., 0.5, 1. ],
    [1.5, 2., 2.5],
    [3., 3.5, 4. ],
    [4.5, 5., 5.5],
    [6., 6.5, 7.]])

两个 ndarray 对象的水平堆叠。

两个 ndarray 对象的垂直堆叠。

另一种特殊操作是将多维 ndarray 对象展平(Flatten)为一维对象。可以选择是按行展平(C 顺序)还是按列展平(F 顺序):

In [90]: h
Out[90]: array([[ 0, 1, 2],
    [ 3, 4, 5],
    [ 6, 7, 8],
    [ 9, 10, 11],
    [12, 13, 14]])
In [91]: h.flatten()
Out[91]: array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])
In [92]: h.flatten(order='C')
Out[92]: array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])
In [93]: h.flatten(order='F')
Out[93]: array([ 0, 3, 6, 9, 12, 1, 4, 7, 10, 13, 2, 5, 8, 11, 14])
In [94]: for i in h.flat:
    print(i, end=',')
    0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,
In [95]: for i in h.ravel(order='C'):
print(i, end=',')
0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,
In [96]: for i in h.ravel(order='F'): ④
print(i, end=',')
0,3,6,9,12,1,4,7,10,13,2,5,8,11,14,

展平的默认顺序是 C 顺序。

使用 F 顺序展平。

flat 属性提供扁平迭代器(C 顺序)。

ravel() 方法是 flatten() 的替代方案。

布尔数组#

比较和逻辑运算通常以逐元素的方式在 ndarray 对象上工作,与标准Python数据类型相同。默认情况下,评估条件产生一个布尔 ndarray 对象(dtype 为 bool):

In [97]: h
Out[97]: array([[0, 1, 2],
    [3, 4, 5],
    [6, 7, 8],
    [9, 10, 11],
    [12, 13, 14]])
In [98]: h > 8 ①
Out[98]: array([[False, False, False],
    [False, False, False],
    [False, False, False],
    [True, True, True],
    [True, True, True]])
In [99]: h <= 7 ②
Out[99]: array([[True, True, True],
    [True, True, True],
    [True, True, False],
    [False, False, False],
    [False, False, False]])
In [100]: h == 5 ③
Out[100]: array([[False, False, False],
    [False, False, True],
    [False, False, False],
    [False, False, False],
    [False, False, False]])
In [101]: (h == 5).astype(int) ④
Out[101]: array([[0, 0, 0],
[0, 0, 1],
[0, 0, 0],
[0, 0, 0],
[0, 0, 0]])
In [102]: (h > 4) & (h <= 12)
Out[102]: array([[False, False, False],
    [False, False, True],
    [True, True, True],
    [True, True, True],
    [True, False, False]])

值是否大于……?

值是否小于或等于……?

值是否等于……?

将 True 和 False 表示为整数值 0 和 1。

值是否大于……且小于或等于……?

这样的布尔数组可用于索引和数据选择。注意以下操作会扁平化数据:

In [103]: h[h > 8] ①
Out[103]: array([9, 10, 11, 12, 13, 14])

In [104]: h[(h > 4) & (h <= 12)] ②
Out[104]: array([5, 6, 7, 8, 9, 10, 11, 12])

In [105]: h[(h < 4) | (h >= 12)] ③
Out[105]: array([0, 1, 2, 3, 12, 13, 14])

给我所有大于……的值。

给我所有大于……且小于或等于……的值。

给我所有大于……或小于或等于……的值。

在这方面一个强大的工具是 np.where() 函数,它允许根据条件为 True 还是 False 来定义操作。应用 np.where() 的结果是一个与原始对象形状相同的新 ndarray 对象:

In [106]: np.where(h > 7, 1, 0)
Out[106]: array([[0, 0, 0],
    [0, 0, 0],
    [0, 0, 1],
    [1, 1, 1],
    [1, 1, 1]])
In [107]: np.where(h % 2 == 0, 'even', 'odd') ②
Out[107]: array([['even', 'odd', 'even'], ['odd', 'even', 'odd'], ['even', 'odd', 'even'], ['odd', 'even', 'odd'], ['even', 'odd', 'even']], dtype='<U4')
In [108]: np.where(h <= 7, h * 2, h / 2) ③
Out[108]: array([0., 2., 4.], [6., 8., 10.], [12., 14., 4.], [4.5, 5., 5.5], [6., 6.5, 7.]])

在新对象中,如果为 True 则设为 1,否则设为 0。

在新对象中,如果为 True 则设为 'even',否则设为 'odd'。

在新对象中,如果为 True 则设为 h 元素的两倍,否则设为 h 元素的一半。

后续章节将提供更多关于这些重要 ndarray 操作的示例。

速度比较#

我们很快会进入 NumPy 的结构化数组部分,但先让我们再关注一下常规数组,看看专门化带来了什么性能优势。

作为一个简单示例,考虑生成一个形状为 5,000 × 5,000 的矩阵/数组,填充伪随机标准正态分布数。然后计算所有元素的和。首先是纯 Python 方法,使用列表推导式:

In [109]: import random
    I = 5000

In [110]: %time mat = [[random.gauss(0, 1) for j in range(I)] \
    for i in range(I)]
    CPU times: user 17.1 s, sys: 361 ms, total: 17.4 s
    Wall time: 17.4 s

In [111]: mat[0][:5] ②
Out[111]: [-0.40594967782329183,
    -1.357757478015285,
    0.05129566894355976,
    -0.8958429976582192,
    0.6234174778878331]
In [114]: %time mat = np.random.standard_normal((I, I)) ①
CPU times: user 1.01 s, sys: 200 ms, total: 1.21 s
Wall time: 1.21 s

In [115]: %time mat.sum() ②
CPU times: user 29.7 ms, sys: 1.15 ms, total: 30.8 ms
Wall time: 29.4 ms

Out[115]: -186.12767026606448

In [116]: mat.nbytes ③
Out[116]: 200000000

In [117]: sys.getsizeof(mat) ③
Out[117]: 200000112
In [112]: %time sum([sum(l) for l in mat]) ③
CPU times: user 142 ms, sys: 1.69 ms, total: 144 ms
Wall time: 143 ms

Out[112]: -3561.944965714259

In [113]: import sys
sum([sys.getsizeof(l) for l in mat]) ④
Out[113]: 215200000

通过嵌套列表推导式创建矩阵。

一些选中的随机数。

在列表推导式中先计算每个列表对象的和,然后对所有和求和。

这累加所有列表对象的内存使用。

现在让我们转向 NumPy,看看相同问题是如何解决的。为方便起见,NumPy 的子包 random 提供了大量函数,可以同时实例化 ndarray 对象并用伪随机数填充它:

创建包含标准正态分布随机数的 ndarray 对象;速度快了约14倍。

计算 ndarray 对象中所有值的和;速度快了4.5倍。

NumPy 方法还节省了一些内存,因为 ndarray 对象的内存开销与数据本身的大小相比微乎其微。

对于基于数组的操作和算法,使用 NumPy 通常会产生紧凑、易读的代码,并比纯 Python 代码有显著的性能提升。

结构化NumPy数组#

ndarray 类的专门化显然带来了许多有价值的优势。然而,过于狭窄的专门化可能会成为大多数基于数组的算法和应用程序的沉重负担。因此,NumPy 提供了结构化 ndarray 和记录 recarray 对象,允许每列有不同的 dtype。"每列"是什么意思?考虑以下结构化 ndarray 对象的初始化:

In [118]: dt = np.dtype([('Name', 'S10'), ('Age', 'i4'),
    ('Height', 'f'), ('Children/Pets', 'i4', 2)]) ①
In [119]: dt ①
Out[119]: dtype([('Name', 'S10'), ('Age', '<i4'), ('Height', '<f4'),
    ('Children/Pets', '<i4', (2,))])
In [120]: dt = np.dtype({'names': ['Name', 'Age', 'Height', 'Children/Pets'],
    'formats': '0 int float int,int'.split()}) ②
In [121]: dt ②
Out[121]: dtype([('Name', 'O'), ('Age', '<i8'), ('Height', '<f8'),
    ('Children/Pets', [('f0', '<i8'), ('f1', '<i8')])])
In [122]: s = np.array([('Smith', 45, 1.83, (0, 1)),
    ('Jones', 53, 1.72, (2, 2))], dtype=dt) ③
In [123]: s ③
Out[123]: array([('Smith', 45, 1.83, (0, 1)), ('Jones', 53, 1.72, (2, 2))],
    dtype=[('Name', 'O'), ('Age', '<i8'), ('Height', '<f8'),
    ('Children/Pets', [('f0', '<i8'), ('f1', '<i8')])])
In [124]: type(s) ④
Out[124]: numpy.ndarray

组合出复杂的 dtype。

另一种语法,达到相同结果。

用两条记录实例化结构化 ndarray。

对象类型仍然是 ndarray。

In [125]: s['Name'] ①
Out[125]: array(['Smith', 'Jones'], dtype=object)
In [126]: s['Height'].mean() ②
Out[126]: 1.775
In [127]: s[0] ③
Out[127]: ('Smith', 45, 1.83, (0, 1))
In [128]: s[1]['Age'] ④
Out[128]: 53

在某种意义上,这种构造非常接近于在SQL数据库中初始化表的操作:有列名和列数据类型,可能还有一些附加信息(例如每个 str 对象的最大字符数)。现在可以按名称轻松访问各列,按索引值访问各行:

按名称选择列。

在选定列上调用方法。

选择一条记录。

选择记录中的一个字段。

总之,结构化数组是常规 ndarray 对象类型的泛化,因为数据类型只需在每列内相同,就像SQL数据库中的表一样。结构化数组的一个优点是列中的单个元素可以是另一个多维对象,不必符合基本的 NumPy 数据类型。

结构化数组#

除了常规数组,NumPy 还提供结构化(和记录)数组,允许描述和处理类似于表的数据结构,每列(有名称)可以有多种不同的数据类型。它们将类似SQL表的数据结构带入Python,同时具有常规 ndarray 对象的大部分优势(语法、方法、性能)。

代码向量化#

向量化(Vectorization)是一种获得更紧凑代码、且可能执行更快的策略。其基本思想是"一次性"对复杂对象进行操作或将函数应用于复杂对象,而不是通过循环遍历该对象的单个元素。在

In [129]: np.random.seed(100)
    r = np.arange(12).reshape((4, 3)) ①
    s = np.arange(12).reshape((4, 3)) * 0.5 ②

In [130]: r ①
Out[130]: array([[0, 1, 2],
    [3, 4, 5],
    [6, 7, 8],
    [9, 10, 11]])
In [131]: s ②
Out[131]: array([[0., 0.5, 1.],
    [1.5, 2., 2.5],
    [3., 3.5, 4.],
    [4.5, 5., 5.5]])
In [132]: r + s ③
Out[132]: array([[0., 1.5, 3.],
    [4.5, 6., 7.5],
    [9., 10.5, 12.],
    [13.5, 15., 16.5]])

Python中,像 map() 和 filter() 这样的函数式编程工具提供了一些基本的向量化手段。然而,NumPy 在其核心深处内置了向量化功能。

基本向量化#

如前一小节所示,简单的数学运算——如计算所有元素的和——可以直接在 ndarray 对象上实现(通过方法或通用函数)。更一般的向量化操作也是可能的。例如,可以按如下方式逐元素相加两个 NumPy 数组:

第一个 ndarray 对象包含随机数。

第二个 ndarray 对象包含随机数。

逐元素相加作为向量化操作(无循环)。

NumPy 还支持所谓的广播(Broadcasting)。这允许在单个操作中组合不同形状的对象。前面的例子已经用到了这一点。考虑以下示例:

In [133]: r + 3
Out[133]: array([[3, 4, 5],
    [6, 7, 8],
    [9, 10, 11],
    [12, 13, 14]])
In [134]: 2 * r
Out[134]: array([[0, 2, 4],
    [6, 8, 10],
    [12, 14, 16],
    [18, 20, 22]])
In [135]: 2 * r + 3
Out[135]: array([[3, 5, 7],
    [9, 11, 13],
    [15, 17, 19],
    [21, 23, 25]])

在标量加法期间,标量被广播并加到每个元素。

在标量乘法期间,标量也被广播并与每个元素相乘。

这个线性变换结合了两种操作。

这些操作在一定程度上也适用于形状不同的 ndarray 对象:

In [136]: r
Out[136]: array([[0, 1, 2],
    [3, 4, 5],
    [6, 7, 8],
    [9, 10, 11]])
In [137]: r.shape
Out[137]: (4, 3)

In [138]: s = np.arange(0, 12, 4) ①
s ①
Out[138]: array([0, 4, 8])

In [139]: r + s ②
Out[139]: array([[0, 5, 10],
    [3, 8, 13],
    [6, 11, 16],
    [9, 14, 19]])
In [140]: s = np.arange(0, 12, 3) ③
s ③
Out[140]: array([0, 3, 6, 9])

In [141]: r + s ④

ValueErrorTraceback (most recent call last)
<ipython-input-141-1890b26ec965> in <module>()
----> 1 r + s ④
ValueError: operands could not be broadcast together with shapes (4,3) (4,)

In [142]: r.transpose() + s ⑤
Out[142]: array([[0, 6, 12, 18],
    [1, 7, 13, 19],
    [2, 8, 14, 20]])
In [143]: sr = s.reshape(-1, 1) ⑥
sr
Out[143]: array([[0],
    [3],
    [6],
    [9]])
In [144]: sr.shape ⑥
Out[144]: (4, 1)

In [145]: r + s.reshape(-1, 1) ⑥
Out[145]: array([[0, 1, 2],
    [6, 7, 8],
    [12, 13, 14],
    [18, 19, 20]])

一个新的长度为3的一维 ndarray 对象。

r(矩阵)和 s(向量)对象可以直接相加。

另一个长度为4的一维 ndarray 对象。

新的 s(向量)对象的长度现在与 r 对象的第二维长度不同。

转置 r 对象后,可以再次进行向量化加法。

或者,将 s 的形状更改为 (4, 1) 可以使加法生效(但结果不同)。

通常,自定义的Python函数也可以与 ndarray 对象一起使用。如果实现允许,数组可以与 int 或 float 对象一样用于函数。考虑以下函数:

In [146]: def f(x):
    return 3 * x + 5
In [147]: f(0.5)
Out[147]: 6.5
In [148]: f(r)
Out[148]: array([[5, 8, 11], [14, 17, 20], [23, 26, 29], [32, 35, 38]])

一个简单的Python函数,实现对参数 x 的线性变换。

函数 f() 应用于Python float 对象。

相同的函数应用于 ndarray 对象,得到向量化的逐元素函数求值。

NumPy 所做的就是简单地逐元素将函数 f 应用于对象。在这个意义上,使用这种操作并不能避免循环;只是避免了 Python 层面的循环,将循环委托给 NumPy。在 NumPy 层面,由优化代码负责遍历 ndarray 对象,这些代码大部分是用 C 编写的,因此通常比纯 Python 更快。这就解释了使用 NumPy 进行基于数组的用例时性能优势背后的"秘密"。

内存布局#

当使用 np.zeros() 初始化 ndarray 对象时(如"多维" on page 94 所示),提供了一个可选参数用于内存布局。这个参数大致指定了数组中的哪些元素在内存中相邻存储(连续)。处理小数组时,这对数组操作的性能几乎没有什么可衡量的影响。然而,当数组变大,并且取决于要实现的(金融)算法时,情况可能会不同。这时内存布局就变得重要了(参见例如 Eli Bendersky 的文章 "Memory Layout of Multi-Dimensional Arrays")。

为了说明数组内存布局在科学和金融领域的潜在重要性,考虑以下多维 ndarray 对象的构造:

In [149]: x = np.random.standard_normal((1000000, 5))
In [150]: y = 2 * x + 3
In [151]: C = np.array((x, y), order='C')
In [152]: F = np.array((x, y), order='F')
In [153]: x = 0.0; y = 0.0
In [154]: C[:2].round(2)
Out[154]: array([[[-1.75, 0.34, 1.15, -0.25, 0.98], [0.51, 0.22, -1.07, -0.19, 0.26],
[-0.46, 0.44, -0.58, 0.82, 0.67],
...
[-0.05, 0.14, 0.17, 0.33, 1.39],
[1.02, 0.3, -1.23, -0.68, -0.87],
[0.83, -0.73, 1.03, 0.34, -0.46]],
[[-0.5, 3.69, 5.31, 2.5, 4.96],
[4.03, 3.44, 0.86, 2.62, 3.51],
[2.08, 3.87, 1.83, 4.63, 4.35],
...
[2.9, 3.28, 3.33, 3.67, 5.78],
[5.04, 3.6, 0.54, 1.65, 1.26],
[4.67, 1.54, 5.06, 3.69, 2.07]]])

一个在两个维度上具有较大不对称性的 ndarray 对象。

原始对象数据的线性变换。

创建一个 C 顺序(行优先)的二维 ndarray 对象。

创建一个 F 顺序(列优先)的二维 ndarray 对象。

释放内存(取决于垃圾回收)。

来自 C 对象的一些数字。

让我们看看这两种 ndarray 对象的一些基本示例和用例,并考虑不同内存布局下的执行速度:

In [155]: %timeit C.sum() ①
4.36 ms ± 89.3 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [156]: %timeit F.sum() ①
4.21 ms ± 71.4 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [157]: %timeit C.sum(axis=0) ②
17.9 ms ± 776 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [158]: %timeit C.sum(axis=1) ③
35.1 ms ± 999 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [159]: %timeit F.sum(axis=0) ②
83.8 ms ± 2.63 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [160]: %timeit F.sum(axis=1) ③
67.9 ms ± 5.16 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [161]: F = 0.0; C = 0.0

计算所有元素的和。

计算每行的和("多个")。

计算每列的和("少量")。

我们可以将性能结果总结如下:

• 计算所有元素的和时,内存布局并不重要。

• C 顺序的 ndarray 对象在行和列上的求和都更快(绝对速度优势)。

• 对于 C 顺序(行优先)的 ndarray 对象,按行求和相对比按列求和更快。

• 对于 F 顺序(列优先)的 ndarray 对象,按列求和相对比按行求和更快。

结语#

NumPy 是 Python 中进行数值计算的首选包。ndarray 类专门设计用于方便高效地处理(大型)数值数据。强大的方法和 NumPy 通用函数允许编写向量化代码,从而避免在 Python 层面上使用缓慢的循环。本章介绍的许多方法也适用于 pandas 及其 DataFrame 类(参见第5章)。

延伸资源#

NumPy 网站提供了许多有用的资源:

• http://www.numpy.org/

推荐的 NumPy 书籍参考资料包括:

• McKinney, Wes (2017). Python for Data Analysis. Sebastopol, CA: O'Reilly.

• VanderPlas, Jake (2016). Python Data Science Handbook. Sebastopol, CA: O'Reilly.