NumPy学习笔记—(1/3)

可以叫我才哥

发布于 2021-08-05 11:21:26

1.5K0

发布于 2021-08-05 11:21:26

文章被收录于专栏：可以叫我才哥

NumPy（Numerical Python）是Python的一种开源的数值计算扩展。这种工具可用来存储和处理大型矩阵，比Python自身的嵌套列表（nested list structure)结构要高效的多（该结构也可以用来表示矩阵（matrix）），支持大量的维度数组与矩阵运算，此外也针对数组运算提供大量的数学函数库。

1.理解 Python 中的数据类型

想要有效的掌握数据驱动科学和计算需要理解数据是如何存储和处理的。本节将描述和对比数组在 Python 语言中和在 NumPy 中是怎么处理的，NumPy 是如何优化了这部分的内容。

Python 的用户通常都是被它的易用性吸引来的，其中很重要一环就是动态类型。静态类型的语言，例如 C 或者 Java，每个变量都需要明确声明，而动态类型语言如 Python 就略过了这个部分。例如，在 C 中，你可能会写如下的代码片段：

int result = 0;
for(int i=0; i<100; i++){
    result += i;
}

但是在 Python 当中，等效的代码如下：

result = 0
for i in range(100):
    result += i

注意其中主要的区别：在 C 当中，每个变量都需要显式声明，Python 的类型是动态推断的。这意味着，我们可以给任何的变量赋值为任何类型的数据，例如：

#python环境下
x = 4
x = "four"

上面的例子中我们将x变量的内容从一个整数变成了一个字符串。如果你想在 C 语言中这样做，取决于不同的编译器，可能会导致一个编译错误或者其他无法预料的结果。

#C语言环境下
int x = 4;
x = "four";  // 编译错误

这种灵活性提供了 Python 和其他动态类型语言在使用上的简易性。但是，理解这里面的工作原理对于在 Python 中高效准确的学习和分析数据是非常重要的。Python 的这种类型灵活性，实际上是付出了额外的存储代价的，变量不仅仅存储了数据本身，还需要存储其相应的类型。我们会在本节接下来的部分继续讨论。

1.1.Python 的整数不仅仅是一个整数

标准的 Python 实现是使用 C 语言编写的。这意味着每个 Python 当中的对象都是一个伪装良好的 C 结构体，结构体内不仅仅包括它的值，还有其他的信息。例如，当我们在 Python 中定义了一个整数，比方说x=10000，x不仅仅是一个原始的整数，它在底层实际上是一个指向复杂 C 结构体的指针，里面含有若干个字段。当你查阅 Python 3.4 的源代码的时候，你会发现整数（实际上是长整形）的定义如下（我们将 C 语言中的宏定义展开后）：

struct _longobject {
    long ob_refcnt;
    PyTypeObject *ob_type;
    size_t ob_size;
    long ob_digit[1];
};

一个 Python 的整数实际上包含四个部分：

ob_refcnt：引用计数器，Python 用这个字段来进行内存分配和垃圾收集
ob_type：变量类型的编码内容
ob_size：表示下面的数据字段的长度
ob_digit：真正的整数值存储在这个字段

这意味着在 Python 中存储一个整数要比在像 C 这样的编译语言中存储一个整数要有损耗，就像下图展示的那样：

这里的PyObject_HEAD代表了前面的引用计数器、类型代码和数据长度的三个字段内容。

再次注意一下这里的区别：C 的整数就是简单一个内存位置，这个位置上的固定长度的字节可以表示一个整数；Python 中的一个整数是一个指向内存位置的指针，该内存位置包括 Python 需要表示一个整数的所有信息，其中最后固定长度的字节才真正存储这个整数。这些额外的信息提供了 Python 的灵活性和易用性。这些 Python 类型需要的额外信息是有额外损失的，特别是当有一个集合需要存储许多这种类型的数据时。

1.2.Python 的列表不仅仅是一个列表

现在我们继续考虑当我们使用 Python 的数据结构来存储许多这样的 Python 对象时的情况。Python 中标准的可变多元素的容器集合就是列表。我们按如下的方式创建一个整数的列表：

L = list(range(10))
L

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

type(L[0])

int

又或者，类似的，字符串的列表：

L2 = [str(c) for c in L] # 列表解析
L2

['0', '1', '2', '3', '4', '5', '6', '7', '8', '9']

type(L2[0])

str

因为 Python 是动态类型，我们甚至可以创建不同类型元素的列表：

L3 = [True, "2", 3.0, 4]
[type(item) for item in L3]

[bool, str, float, int]

这种灵活性是要付出代价的：要让列表能够容纳不同的类型，每个列表中的元素都必须带有自己的类型信息、引用计数器和其他的信息，一句话，里面的每个元素都是一个完整的 Python 的对象。如果在所有的元素都是同一种类型的情况下，这里面绝大部分的信息都是冗余的：如果我们能将数据存储在一个固定类型的数组中，显然会更加高效。下图展示了动态类型的列表和固定类型的数组（NumPy 实现的）的区别：

Array Memory Layout

从底层实现上看，数组仅仅包含一个指针指向一块连续的内存空间。而 Python 列表，含有一个指针指向一块连续的指针内存空间，里面的每个指针再指向内存中每个独立的 Python 对象，如我们前面看到的整数。列表的优势在于灵活：因为每个元素都是完整的 Python 的类型对象结构，包含了数据和类型信息，因此列表可以存储任何类型的数据。NumPy 使用的固定类型的数组缺少这种灵活性，但是对于存储和操作数据会高效许多。

1.3.Python 的固定类型数组

Python 提供了许多不同的选择能让你高效的存储数据，使用固定类型数据。內建的array模块（从 Python 3.3 开始提供）可以用来创建同一类型的数组：

import array
L = list(range(10))
A = array.array('i', L)
A

array('i', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

这里的i是表示数据类型是整数的类型代码。

更常用的是ndarray对象，由 NumPy 包提供。虽然 Python 的array提供了数组的高效存储，NumPy 则提供了数组的更高效运算。我们会在后续小节中陆续介绍这些操作，这里我们首先介绍创建 NumPy 数组的方式。

当然最开始要做的是将 NumPy 包载入，惯例上提供别名np：

import numpy as np

1.4.使用 Python 列表创建数组

首先，我们可以使用np.array来将一个 Python 列表变成一个数组：

# 整数数组:
np.array([1, 4, 2, 5, 3])

array([1, 4, 2, 5, 3])

记住和 Python 列表不同，NumPy 数组只能含有同一种类型的数据。如果类型不一样，NumPy 会尝试向上扩展类型（下面例子中会将整数向上扩展为浮点数）：

np.array([3.14, 4, 2, 3])

array([3.14, 4.  , 2.  , 3.  ])

如果你需要明确指定数据的类型，你可以使用dtype关键字参数：

np.array([1, 2, 3, 4], dtype='float32')

array([1., 2., 3., 4.], dtype=float32)

最后，不同于 Python 的列表，NumPy 的数组可以明确表示为多维；下面例子是一个使用列表的列表来创建二维数组的方法：

# 更准确的说，应该是生成器的列表，列表解析中有三个range生成器
# 分别是range(2, 5), range(4, 7) 和 range(6, 9)
np.array([range(i, i + 3) for i in [2, 4, 6]])

array([[2, 3, 4],
       [4, 5, 6],
       [6, 7, 8]])

内部的列表作为二维数组的行。

1.5.从头开始创建数组

使用 NumPy 的方法从头创建数组会更加高效，特别对于大型数组来说。下面有几个例子：

# zeros将数组元素都填充为0，10是数组长度
np.zeros(10, dtype=int)

array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

# ones将数组元素都填充为1，(3, 5)是数组的维度说明，表明数组是二维的3行5列
np.ones((3, 5), dtype=float)

array([[1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.]])

# full将数组元素都填充为参数值3.14，(3, 5)是数组的维度说明，表明数组是二维的3行5列
np.full((3, 5), 3.14)

array([[3.14, 3.14, 3.14, 3.14, 3.14],
       [3.14, 3.14, 3.14, 3.14, 3.14],
       [3.14, 3.14, 3.14, 3.14, 3.14]])

# arange类似range，创建一段序列值
# 起始值是0（包含），结束值是20（不包含），步长为2
np.arange(0, 20, 2)

array([ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18])

# linspace创建一段序列值，其中元素按照区域进行线性（平均）划分
# 起始值是0（包含），结束值是1（包含），共5个元素
np.linspace(0, 1, 5)

array([0.  , 0.25, 0.5 , 0.75, 1.  ])

# random.random随机分布创建数组
# 随机值范围为[0, 1)，(3, 3)是维度说明，二维数组3行3列
np.random.random((3, 3))

array([[0.83591315, 0.97973367, 0.87453881],
       [0.28872419, 0.49923571, 0.66769524],
       [0.51712704, 0.16038453, 0.20525946]])

# random.normal正态分布创建数组
# 均值0，标准差1，(3, 3)是维度说明，二维数组3行3列
np.random.normal(0, 1, (3, 3))

array([[ 1.66414086, -1.35090737,  0.54766935],
       [ 1.43521621, -0.79276737, -0.97960272],
       [-0.00446251,  1.33370298,  0.31689258]])

# random.randint随机整数创建数组，随机数范围[0, 10)
np.random.randint(0, 10, (3, 3))

array([[2, 9, 2],
       [6, 1, 9],
       [2, 4, 8]])

# 3x3的单位矩阵数组
np.eye(3)

array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

# empty创建一个未初始化的数组，数组元素的值保持为原有的内存空间值
np.empty(3)

array([1., 1., 1.])

1.6.NumPy 标准数据类型

NumPy 数组仅包含一种类型数据，因此它的类型系统和 Python 也有所区别，因为对于每一种 NumPy 类型，都需要更详细的类型信息和限制。因为 NumPy 是使用 C 构建的，它的类型系统对于 C、Fortran 的用户来说不会陌生。

标准 NumPy 数据类型见下表。正如上面介绍的，当我们创建数组的时候，我们可以将dtype参数指定为下面类型的字符串名称来指定数组的数据类型。

np.zeros(10, dtype='int16')

也可以将dtype指定为对应的 NumPy 对象：

np.zeros(10, dtype=np.int16)

Data type	Description
bool_	布尔(True 或 False) 一个字节
int_	默认整数类型 (类似 C 的long; 通常可以是int64或int32)
intc	类似 C 的int (通常可以是int32或int64)
intp	用于索引值的整数(类似 C 的ssize_t; 通常可以是int32或int64)
int8	整数，1 字节 (-128 ～ 127)
int16	整数，2 字节 (-32768 ～ 32767)
int32	整数，4 字节 (-2147483648 ～ 2147483647)
int64	整数，8 字节 (-9223372036854775808 ～ 9223372036854775807)
uint8	字节 (0 ～ 255)
uint16	无符号整数 (0 ～ 65535)
uint32	无符号整数 (0 ～ 4294967295)
uint64	无符号整数 (0 ～ 18446744073709551615)
float_	float64的简写
float16	半精度浮点数: 1 比特符号位, 5 比特指数位, 10 比特尾数位
float32	单精度浮点数: 1 比特符号位, 8 比特指数位, 23 比特尾数位
float64	双精度浮点数: 1 比特符号位, 11 比特指数位, 52 比特尾数位
complex_	complex128的简写
complex64	复数, 由 2 个单精度浮点数组成
complex128	复数, 由 2 个双精度浮点数组成

还有更多的高级的类型声明，比如指定大尾或小尾表示；需要获得更多内容，请查阅NumPy 在线文档[1]。NumPy 也支持复合数据类型，这部分我们将在结构化数据：NumPy 里的结构化数组中进行介绍

2.NumPy 数组基础

Python 中的数据操作基本就是 NumPy 数组操作的同义词：一些新的工具像 Pandas 都是依赖于 NumPy 数组建立起来的。本节会展示使用 NumPy 数组操作和访问数据以及子数组的一些例子，包括切分、变形和组合。尽管这里展示的操作有些枯燥和学术化，但是它们是组成本书后面使用的例子的基础。你应该更好的掌握它们。

我们会讨论下述数组操作的基本内容：

数组的属性: 获得数组的大小、形状、内存占用以及数据类型
数组索引: 获得和设置单个数组元素的值
数组切片: 获得和设置数组中的子数组
数组变形: 改变数组的形状
组合和切分数组: 将多个数组组合成一个，或者将一个数组切分成多个

2.1.NumPy 数组属性

首先我们来讨论一些数组有用的属性。我们从定义三个数组开始，一个一维的，一个二维的和一个三维的数组。我们采用 NumPy 的随机数产生器来创建数组，产生之前我们会给定一个随机种子，这样来保证每次代码运行的时候都能得到相同的数组：

import numpy as np
np.random.seed(0)  # 设定随机种子，保证实验的可重现

x1 = np.random.randint(10, size=6)  # 一维数组
x2 = np.random.randint(10, size=(3, 4))  # 二维数组
x3 = np.random.randint(10, size=(3, 4, 5))  # 三维数组

每个数组都有属性ndim，代表数组的维度，shape代表每个维度的长度（形状）和size代表数组的总长度（元素个数）

# 输出三维数组的维度、形状和总长度
print("x3 ndim: ", x3.ndim)
print("x3 shape:", x3.shape)
print("x3 size: ", x3.size)

x3 ndim:  3
x3 shape: (3, 4, 5)
x3 size:  60

另一个有用的属性是dtype，数组的数据类型（我们在上一节理解 Python 的数据类型中已经见过）。

print("dtype:", x3.dtype)

dtype: int32

还有属性包括itemsize代表每个数组元素的长度（单位字节），nbytes代表数组的总字节长度：

print(f'itemsize: {x3.itemsize} bytes')#f-string输出
print("nbytes:", x3.nbytes, "bytes")

itemsize: 4 bytes
nbytes: 240 bytes

通常，我们可以认为nbytes等于itemsize乘以size。

2.2.数组索引：获取单个元素

如果我们熟悉 Python 列表的索引方式，那么 NumPy 数组的索引方式也是很相似的。对于一维数组来说，第 i 个元素值（从 0 开始）可以使用中括号内的索引值获得：

x1

array([5, 0, 3, 3, 7, 9])

x1[0]

x1[4]

需要从末尾进行索引取值，你可以使用负的索引值：

x1[-1]

x1[-2]

在多维数组中获取元素值，可以在中括号中使用一个索引值的元组：

多维数组的索引方式与列表的列表索引方式是不同的。列表的列表在 Python 中需要使用多个中括号进行索引，如x[i][j]的方式。

x2

array([[3, 5, 2, 4],
       [7, 6, 8, 8],
       [1, 6, 7, 7]])

x2[0, 0]

x2[2]

array([1, 6, 7, 7])

x2[2][0]

x2[2, 0]

x2[2, -1]

元素值也可以通过上述的索引语法进行修改：

x2[0, 0] = 12
x2

array([[12,  5,  2,  4],
       [ 7,  6,  8,  8],
       [ 1,  6,  7,  7]])

请记住，与 Python 的列表不同，NumPy 数组是固定类型的。这意味着，如果你试图将一个浮点数值放入一个整数型数组，这个值会被默默地截成整数。这是比较容易犯的错误。

x1[0] = 3.14159  # 会被截成整数
x1

array([3, 0, 3, 3, 7, 9])

2.3.数组切片：获取子数组

x[start:stop:step]

正如我们可以使用中括号获取单个元素值，我们也可以使用中括号的切片语法获取子数组，切片的语法遵从标准 Python 列表的切片语法格式；对于一个数组x进行切片：

x[start:stop:step]

如果三个参数没有设置值的话，默认值分别是start=0，stop=维度的长度，step=1。我们来看看在一维数组和多维数组中进行切片取子数组的例子。

2.3.1.一维子数组

x = np.arange(10)
x

array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

x[:5]  # 前五个元素

array([0, 1, 2, 3, 4])

x[5:]  # 从序号5开始的所有元素

array([5, 6, 7, 8, 9])

x[4:7]  # 中间4~6序号的元素

array([4, 5, 6])

x[::2]  # 每隔一个取元素

array([0, 2, 4, 6, 8])

x[1::2]  # 每隔一个取元素，开始序号为1

array([1, 3, 5, 7, 9])

当 step 为负值时，将会在数组里反向的取元素，这是将数组反向排序最简单的方法：

从其他编程语言转 Python 的初学者，很容易问一个问题，我想反序一个字符串，怎么找不到函数啊，內建的没有，str 的方法也没有。答案是，因为根本不需要，例如： ”

s = 'hello world'
# 下面就会输出'dlrow olleh'
print(s[::-1])

x[::-1]  # 反序数组

array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])

x[5::-2]  # 从序号5开始向前取元素，每隔一个取一个元素

array([5, 3, 1])

2.3.2.多维子数组

多维数组的切片也一样，只是在中括号中使用逗号分隔多个切片声明。例如：

x2

array([[12,  5,  2,  4],
       [ 7,  6,  8,  8],
       [ 1,  6,  7,  7]])

x2[:2, :3]  # 行的维度取前两个，列的维度取前三个，形状变为(2, 3)

array([[12,  5,  2],
       [ 7,  6,  8]])

x2[:3, ::2]  # 行的维度取前三个（全部），列的维度每个一个取一列，形状变为(3, 2)

array([[12,  2],
       [ 7,  8],
       [ 1,  7]])

最后，子数组的各维度还可以反序：

x2[::-1, ::-1] # 行和列都反序，形状保持(3, 4)

array([[ 7,  7,  6,  1],
       [ 8,  8,  6,  7],
       [ 4,  2,  5, 12]])

2.3.2.1.获取数组的行和列

还有一种常见的需要是获取数组的单行或者单列。这可以通过组合索引和切片两个操作做到，使用一个不带参数的冒号:可以表示取该维度的所有元素：

print(x2[:, 0])  # x2的第一列

[12  7  1]

print(x2[0, :])  # x2的第一行

[12  5  2  4]

如果是获取行数据的话，可以省略后续的切片，写成更加简洁的方式：

print(x2[0])  # 等同于 x2[0, :]

[12  5  2  4]

2.3.3.子数组是非副本的视图

一个非常重要和有用的概念你需要知道的就是数组的切片返回的实际上是子数组的视图而不是它们的副本。这是 NumPy 数组的切片和 Python 列表的切片的主要区别，列表的切片返回的是副本。用上面的二维做例子：

print(x2)

[[12  5  2  4]
 [ 7  6  8  8]
 [ 1  6  7  7]]

让我们从中取一个

2 \times 2

的子数组：

x2_sub = x2[:2, :2]
print(x2_sub)

[[12  5]
 [ 7  6]]

如果我们修改这个子数组，我们看到原来的数组也会随之更改：

x2_sub[0, 0] = 99
print(x2_sub)

[[99  5]
 [ 7  6]]

print(x2)

[[99  5  2  4]
 [ 7  6  8  8]
 [ 1  6  7  7]]

这个默认行为是很有用的：这意味着当我们在处理大数据集时，我们可以获取和处理其中的部分子数据集而不需要在内存中复制一份数据的副本。

2.3.4.创建数组的副本

尽管使用视图有上述的优点，有时候我们还是需要从数组中复制一份子数组出来。这可以使用copy方法简单的办到：

x2_sub_copy = x2[:2, :2].copy()
print(x2_sub_copy)

[[99  5]
 [ 7  6]]

现在如果我们改变这个子数组，原数组会保持不变：

x2_sub_copy[0, 0] = 42
print(x2_sub_copy)

[[42  5]
 [ 7  6]]

print(x2)

[[99  5  2  4]
 [ 7  6  8  8]
 [ 1  6  7  7]]

2.4.改变数组的形状

另一个数组的常用操作是改变形状。最方便的方式是使用reshape方法实现。例如，如果你希望将 1~9 的数放入一个

3 \times 3

的数组里面，你可以这样做：

grid = np.arange(1, 10).reshape((3, 3))
print(grid)

[[1 2 3]
 [4 5 6]
 [7 8 9]]

注意，改变形状要能成功，原始数组和新的形状的数组的总长度size必须一样。当可能的情况下，reshape会尽量使用原始数组的视图，但是如果原始数组的数据存储在不连续的内存区，就会进行复制。

另外一个常用的改变形状的操作就是将一个一维数组变成二维数组中的一行或者一列。这也可以使用reshape方法实现，或者更简单的方式是使用切片语法中的newaxis属性增加一个维度：

x = np.array([1, 2, 3])

# 使用reshape变为 (1, 3)
x.reshape((1, 3))

array([[1, 2, 3]])

# 使用newaxis，增加行维度，形状也是 (1, 3)
x[np.newaxis, :]

array([[1, 2, 3]])

# 使用reshape变为 (3, 1)
x.reshape((3, 1))

array([[1],
       [2],
       [3]])

# 使用newaxis增加列维度，形状也是 (3, 1)
x[:, np.newaxis]

array([[1],
       [2],
       [3]])

2.5.数组的连接和切分

前面的方法都是在单个数组上进行操作。我们也可以将多个数组组成一个，或者反过来将一个数组切分成多个。下面我们来看看这些操作。

2.5.1 连接数组

在 NumPy 中连接或者组合多个数组，有三个不同的方法np.concatenate，np.vstack和np.hstack。np.concatenate接受一个数组的元组或列表作为第一个参数，如下：

x = np.array([1, 2, 3])
y = np.array([3, 2, 1])
np.concatenate([x, y])

array([1, 2, 3, 3, 2, 1])

你也可以一次连接两个以上的数组：

z = [99, 99, 99]
print(np.concatenate([x, y, z]))

[ 1  2  3  3  2  1 99 99 99]

也可以用来连接二维数组：

grid = np.array([[1, 2, 3],
                 [4, 5, 6]])

# 沿着第一个维度进行连接，即按照行连接，axis=0
np.concatenate([grid, grid])

array([[1, 2, 3],
       [4, 5, 6],
       [1, 2, 3],
       [4, 5, 6]])

# 沿着第二个维度进行连接，即按照列连接，axis=1
np.concatenate([grid, grid], axis=1)

array([[1, 2, 3, 1, 2, 3],
       [4, 5, 6, 4, 5, 6]])

进行连接的数组如果具有不同的维度，使用np.vstack（垂直堆叠）和np.hstack（水平堆叠）会更加清晰：

x = np.array([1, 2, 3])
grid = np.array([[9, 8, 7],
                 [6, 5, 4]])

# 沿着垂直方向进行堆叠
np.vstack([x, grid])

array([[1, 2, 3],
       [9, 8, 7],
       [6, 5, 4]])

# 沿着水平方向进行堆叠
y = np.array([[99],
              [99]])
np.hstack([grid, y])

array([[ 9,  8,  7, 99],
       [ 6,  5,  4, 99]])

类似的，np.dstack会沿着第三个维度（深度）进行堆叠。

2.5.2 切分数组

连接的反操作是切分，主要的方法包括np.split，np.hsplit和np.vsplit。我们可以传递一个列表参数表示切分的点：

x = [1, 2, 3, 99, 99, 3, 2, 1]
x1, x2, x3 = np.split(x, [3, 5]) # 在序号3和序号5处进行切分，返回三个数组
print(x1, x2, x3)

[1 2 3] [99 99] [3 2 1]

你应该已经发现 N 个切分点会返回 N+1 个子数组。相应的np.hsplit和np.vsplit也是相似的：

grid = np.arange(16).reshape((4, 4))
grid

array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

upper, lower = np.vsplit(grid, [2]) # 沿垂直方向切分，切分点行序号为2
print(upper)
print(lower)

[[0 1 2 3]
 [4 5 6 7]]
[[ 8  9 10 11]
 [12 13 14 15]]

left, right = np.hsplit(grid, [2]) # 沿水平方向切分数组，切分点列序号为2
print(left)
print(right)

[[ 0  1]
 [ 4  5]
 [ 8  9]
 [12 13]]
[[ 2  3]
 [ 6  7]
 [10 11]
 [14 15]]

同样np.dsplit会沿着第三个维度切分数组。

3.NumPy 数组运算：通用函数

直到目前为止，我们已经讨论了一些 NumPy 的基本构件；在下面几个小节中，我们会深入讨论 NumPy 能在 Python 数据科学中占据重要地位的原因。简而言之，NumPy 提供了简单和灵活的接口来对数组数据计算进行优化。

对 NumPy 的数组进行计算相较其他普通的实现方式而言是非常快的。快的原因关键在于使用了向量化的操作，因为它们都是通过 NumPy 的通用函数（ufuncs）实现的。希望通过本节的介绍，能让读者习惯使用 ufuncs，它们能使在数组元素上的重复计算更加快速和高效。本节还会介绍许多 NumPy 中最常用的 ufuncs 数学计算方法。

3.1.循环，慢的实现

Python 的默认实现（被称为 CPython）对于一些操作执行效率很低。这部分归咎于语言本身的动态和解释执行特性：因为类型是动态的，因此不到执行时，无法预知变量的类型，因此不能像 C 或者 Fortran 那样预先将代码编译成机器代码来执行。近年来，也出现了很多尝试来弥补这个缺陷：其中比较流行和著名的包括PyPy[2]，Python 的 JIT 编译实现；Cython[3]，可以将 Python 代码转换为可编译的 C 代码；和Numba[4]，可以将 Python 代码片段转换为 LLVM 字节码。

Python 另一个表现相对低效的方面是当重复进行很多细微操作时，比方说对一个数组中的每个元素进行循环操作。例如，我们有一个数组，现在我们需要计算每个元素的倒数。一个很直接的实现方式就像下面的代码：

import numpy as np
np.random.seed(0)

def compute_reciprocals(values):
    output = np.empty(len(values))
    for i in range(len(values)):
        output[i] = 1.0 / values[i]
    return output

values = np.random.randint(1, 10, size=5)
compute_reciprocals(values)

array([0.16666667, 1.        , 0.25      , 0.25      , 0.125     ])

上面的代码实现对于很多具有 C 或者 Java 语言背景的读者来说是非常自然的。但是如果我们在一个很大的数据集上测量上面代码的执行时间，我们会发现这个操作很慢，甚至慢的让你吃惊。下面使用%timeit魔术指令对一个大数据集进行测时：

big_array = np.random.randint(1, 100, size=1000000)
%timeit compute_reciprocals(big_array)

3.01 s ± 159 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

这个操作对于百万级的数据集耗时需要几秒。当现在手机的每秒浮点数运算次数都已经已经达到 10 亿级别，这实在是不可思议的慢了。通过分析发现瓶颈并不是代码本身，而是每次循环时 CPython 必须执行的类型检查和函数匹配。每次计算倒数时，Python 首先需要检查对象的类型，然后寻找一个最合适的函数对这种类型进行计算。如果我们使用编译型的语言实现上面的代码，每次计算的时候，类型和应该执行的函数都已经确定，因此执行的时间肯定短很多。

3.2.UFuncs 介绍

对于许多操作，NumPy 都为这种静态类型提供了编译好的函数。被称为向量化的操作。向量化操作可以简单应用在数组上，实际上会应用在每一个元素上。实现原理就是将循环的部分放进 NumPy 编译后的那个层次，从而提高性能。

比较一下下述两种方式得到的结果：

print(compute_reciprocals(values))
print(1.0 / values)

[0.16666667 1.         0.25       0.25       0.125     ]
[0.16666667 1.         0.25       0.25       0.125     ]

下面使用 ufuncs 来测算执行时间，我们可以看到执行时间相差了好几个数量级：

%timeit (1.0 / big_array)

6.96 ms ± 184 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

NumPy 中的向量化操作是通过ufuncs实现的，其主要目的就是在 NumPy 数组中快速执行重复的元素操作。Ufuncs 是极端灵活的，我们上面看到是标量和数组间的操作，但是我们也可以将它们用在两个数组之间：

np.arange(5) / np.arange(1, 6)

array([0.        , 0.5       , 0.66666667, 0.75      , 0.8       ])

而且 ufuncs 也不仅限于一维数组，多维数组同样适用：

x = np.arange(9).reshape((3, 3))
2 ** x

array([[  1,   2,   4],
       [  8,  16,  32],
       [ 64, 128, 256]], dtype=int32)

通过 ufuncs 向量化计算基本上都会比使用 Python 循环实现的相同方法要更加高效，特别是数组的长度增长的情况下。任何情况下，如果你看到 Python 的数组循环操作，都可以替换成为向量化形式。

3.3.NumPy 的 UFuncs

Ufuncs 有两种类型：一元 ufuncs（仅对一个输入值进行操作）和二元 ufuncs（对两个输入值进行操作）。下面我们会看到它们的使用例子。

3.3.1.数组运算

NumPy 的 ufuncs 用起来非常的自然和人性化，因为它们采用了 Python 本身的算术运算符号 - 标准的加减乘除的实现：

x = np.arange(4)
print("x     =", x)
print("x + 5 =", x + 5)
print("x - 5 =", x - 5)
print("x * 2 =", x * 2)
print("x / 2 =", x / 2)
print("x // 2 =", x // 2)  # 整除

x     = [0 1 2 3]
x + 5 = [5 6 7 8]
x - 5 = [-5 -4 -3 -2]
x * 2 = [0 2 4 6]
x / 2 = [0.  0.5 1.  1.5]
x // 2 = [0 0 1 1]

下面是一元的取反，**求幂和%取模：

print("-x     = ", -x)
print("x ** 2 = ", x ** 2)
print("x % 2  = ", x % 2)

-x     =  [ 0 -1 -2 -3]
x ** 2 =  [0 1 4 9]
x % 2  =  [0 1 0 1]

当然，你可以将这些运算按照你的需要组合起来，运算顺序与标准运算一致：

-(0.5*x + 1) ** 2

array([-1.  , -2.25, -4.  , -6.25])

上面看到的这些算术运算操作，都是 NumPy 中相应函数的简化写法；例如+号实际上是add函数的封装：

np.add(x, 2)

array([2, 3, 4, 5])

下表列出 NumPy 实现的运算符号及对应的 ufunc 函数：

运算符	对应的 ufunc 函数	说明
+	np.add	加法 (例如 1 + 1 = 2)
-	np.subtract	减法 (例如 3 - 2 = 1)
-	np.negative	一元取负 (例如 -2)
*	np.multiply	乘法 (例如 2 * 3 = 6)
/	np.divide	除法 (例如 3 / 2 = 1.5)
//	np.floor_divide	整除 (例如 3 // 2 = 1)
**	np.power	求幂 (例如 2 ** 3 = 8)
%	np.mod	模除 (例如 9 % 4 = 1)

除此之外还有布尔和二进制位操作；我们会在[比较，遮盖和布尔逻辑]中介绍它们。

3.3.2.绝对值

就像 NumPy 能够理解 Python 內建的算术操作一样，它同样能理解 Python 內建的绝对值函数：

x = np.array([-2, -1, 0, 1, 2])
abs(x)

array([2, 1, 0, 1, 2])

对应的 NumPy 的 ufunc 是np.absolute，还有一个简短的别名np.abs：

np.absolute(x)

array([2, 1, 0, 1, 2])

np.abs(x)

array([2, 1, 0, 1, 2])

这个 ufunc 可以处理复数，返回的是矢量的长度：

x = np.array([3 - 4j, 4 - 3j, 2 + 0j, 0 + 1j])
np.abs(x)

array([5., 5., 2., 1.])

3.3.3.三角函数

NumPy 提供了大量的有用的 ufuncs，对于数据科学加来说非常有用的还包括三角函数。我们先定义一个角度的数组：

theta = np.linspace(0, np.pi, 3)

然后来计算这个数组的一些三角函数值：

print("theta      = ", theta)
print("sin(theta) = ", np.sin(theta)) # 正弦
print("cos(theta) = ", np.cos(theta)) # 余弦
print("tan(theta) = ", np.tan(theta)) # 正切

theta      =  [0.         1.57079633 3.14159265]
sin(theta) =  [0.0000000e+00 1.0000000e+00 1.2246468e-16]
cos(theta) =  [ 1.000000e+00  6.123234e-17 -1.000000e+00]
tan(theta) =  [ 0.00000000e+00  1.63312394e+16 -1.22464680e-16]

计算得到的值受到计算机浮点数精度的限制，因为上面看到的结果中应该为 0 的地方并不精确的等于 0。还提供了逆三角函数：

x = [-1, 0, 1]
print("x         = ", x)
print("arcsin(x) = ", np.arcsin(x)) # 反正弦
print("arccos(x) = ", np.arccos(x)) # 反余弦
print("arctan(x) = ", np.arctan(x)) # 反正切

x         =  [-1, 0, 1]
arcsin(x) =  [-1.57079633  0.          1.57079633]
arccos(x) =  [3.14159265 1.57079633 0.        ]
arctan(x) =  [-0.78539816  0.          0.78539816]

3.3.4.指数和对数

NumPy 中另一种常用操作是指数：

x = [1, 2, 3]
print("x     =", x)
print("e^x   =", np.exp(x))
print("2^x   =", np.exp2(x))
print("3^x   =", np.power(3, x))

x     = [1, 2, 3]
e^x   = [ 2.71828183  7.3890561  20.08553692]
2^x   = [2. 4. 8.]
3^x   = [ 3  9 27]

指数的逆操作，对数函数。np.log求的是自然对数；如果你需要计算 2 的对数或者 10 的对数，也有相应的函数：

x = [1, 2, 4, 10]
print("x        =", x)
print("ln(x)    =", np.log(x))
print("log2(x)  =", np.log2(x))
print("log10(x) =", np.log10(x))

x        = [1, 2, 4, 10]
ln(x)    = [0.         0.69314718 1.38629436 2.30258509]
log2(x)  = [0.         1.         2.         3.32192809]
log10(x) = [0.         0.30103    0.60205999 1.        ]

还有当输入值很小时，可以保持精度的指数和对数函数：

x = [0, 0.001, 0.01, 0.1]
print("exp(x) - 1 =", np.expm1(x))
print("log(1 + x) =", np.log1p(x))

exp(x) - 1 = [0.         0.0010005  0.01005017 0.10517092]
log(1 + x) = [0.         0.0009995  0.00995033 0.09531018]

当x很小时，这些函数会比np.log或np.exp计算得到更加精确的结果。

3.3.5.特殊的 ufuncs

NumPy 包含更多的 ufuncs，包括双曲函数，二进制位运算，比较操作，角度弧度转换，舍入以及求余数等等。参考 NumPy 的在线文档你可以看到很多有趣的函数说明。

在scipy.special模块中还有更多的特殊及难懂的 ufuncs。如果你需要计算使用到晦涩数学函数操作你的数据，基本上你都可以在这个模块中找到。下面列出了部分与数据统计相关的 ufuncs，还有很多因为篇幅关系并未列出。

from scipy import special

# 伽玛函数（通用阶乘函数）及相关函数
x = [1, 5, 10]
print("gamma(x)     =", special.gamma(x)) # 伽玛函数
print("ln|gamma(x)| =", special.gammaln(x)) # 伽玛函数的自然对数
print("beta(x, 2)   =", special.beta(x, 2)) # 贝塔函数（第一类欧拉积分）

gamma(x)     = [1.0000e+00 2.4000e+01 3.6288e+05]
ln|gamma(x)| = [ 0.          3.17805383 12.80182748]
beta(x, 2)   = [0.5        0.03333333 0.00909091]

# 误差函数 (高斯函数积分)
# 互补误差函数，逆误差函数
x = np.array([0, 0.3, 0.7, 1.0])
print("erf(x)  =", special.erf(x)) # 误差函数
print("erfc(x) =", special.erfc(x)) # 互补误差函数
print("erfinv(x) =", special.erfinv(x)) # 逆误差函数

erf(x)  = [0.         0.32862676 0.67780119 0.84270079]
erfc(x) = [1.         0.67137324 0.32219881 0.15729921]
erfinv(x) = [0.         0.27246271 0.73286908        inf]

还有很多很多 ufuncs，你可以在 NumPy 和scipy.special中找到。因为这些函数的文档都有在线版本，你可以用"gamma 函数 python"就可以找到相关的信息。

3.4.高级 Ufunc 特性

许多 NumPy 用户在使用 ufuncs 的时候都没有了解它们完整特性。我们在这里会简单介绍一些特别的特性。

3.4.1.指定输出

对于大数据量的计算，有时指定存储输出数据的数组是很有用的。指定输出结果的内存位置能够避免创建临时的数组。所有的 ufuncs 都能通过指定out参数来指定输出的数组。

x = np.arange(5)
y = np.empty(5)
np.multiply(x, 10, out=y) # 指定结果存储在y数组中
print(y)

[ 0. 10. 20. 30. 40.]

输出结果甚至可以指定为数组的视图。例如，你可以将结果隔一个元素写入到一个数组中：

y = np.zeros(10)
np.power(2, x, out=y[::2]) # 指定结果存储在y数组中，每隔一个元素存一个
print(y)

[ 1.  0.  2.  0.  4.  0.  8.  0. 16.  0.]

如果你没使用out参数，而是写成y[::2] = 2 ** x，这回导致首先创建一个临时数组用来存储2 ** x，然后再将这些值复制到 y 数组中。对于上面这么小的数组来说，其实没有什么区别，但是如果对象是一个非常大的数组，使用out参数能节省很多内存空间。

3.4.2.聚合

对于二元运算 ufuncs 来说，还有一些很有趣的聚合函数可以直接从数组中计算出结果。例如，如果你想reduce一个数组，你可以对于任何 ufuncs 应用reduce方法。reduce 会重复在数组的每一个元素进行 ufunc 的操作，直到最后得到一个标量。

例如，在add ufunc 上调用reduce会返回所有元素的总和：

x = np.arange(1, 6)
np.add.reduce(x)

相应的，在multiply ufunc 上调用reduce会返回所有元素的乘积：

np.multiply.reduce(x)

如果你需要得到每一步计算得到的中间结果，你可以调用accumulate：

np.add.accumulate(x)

array([ 1,  3,  6, 10, 15], dtype=int32)

np.multiply.accumulate(x)

array([  1,   2,   6,  24, 120], dtype=int32)

注意对于上面这种特殊情况，NumPy 也提供了相应的函数直接计算结果（np.sum，np.prod，np.cumsum，np.cumprod），我们会在[聚合：Min, Max, 以及其他]中详细讨论。

3.4.3.外积

最后，任何 ufunc 都可以计算输入的每一对元素的结果，使用outer方法。你可以一行代码就完成类似创建乘法表的功能：

x = np.arange(1, 6)
np.multiply.outer(x, x)

array([[ 1,  2,  3,  4,  5],
       [ 2,  4,  6,  8, 10],
       [ 3,  6,  9, 12, 15],
       [ 4,  8, 12, 16, 20],
       [ 5, 10, 15, 20, 25]])

ufunc.at和ufunc.reduceat方法也非常有用，我们会在后面高级索引中详细讨论。

Ufuncs 还有一个极端有用的特性，能让 ufuncs 在不同长度和形状的数组之间进行计算，这是一组被称为广播的方法。这是一个非常重要的内容，因此我们会专门在后面*在数组上计算：广播小节中进行介绍。

3.5.Ufuncs：更多资源

更多有关 ufuncs 的信息（包括完整的函数列表）可以在NumPy[5] 和 SciPy[6]的在线文档获得。

不要忘记了我们可以使用 IPython 的帮助工具?来获取任何相关的帮助信息，。

参考资料

[1]

NumPy在线文档: http://numpy.org/

[2]

PyPy: http://pypy.org/

[3]

Cython: http://cython.org

[4]

Numba: http://numba.pydata.org/

[5]

NumPy: http://www.numpy.org

[6]

SciPy: http://www.scipy.org

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-07-18，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

python

numpy

存储

本文分享自可以叫我才哥微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！