9.5 NumPy 数组上的计算:通用函数
本节是《Python 数据科学手册》(Python Data Science Handbook)的摘录。 译者:飞龙 协议:CC BY-NC-SA 4.0
到目前为止,我们一直在讨论 NumPy 的一些基本要点;在接下来的几节中,我们将深入探讨 NumPy 在 Python 数据科学领域如此重要的原因。也就是说,它为数据数组的最优计算,提供了一个简单而灵活的接口。
NumPy 数组的计算速度非常快,也可能非常慢。使其快速的关键是使用向量化操作,通常通过 NumPy 的通用函数(ufunc
)实现。
本节激发了 NumPy 的ufunc
的需求,这些ufunc
可用于更有效地对数组元素进行重复计算。然后介绍了 NumPy 包中可用的,许多最常用和最有用的算术ufunc
。
慢速的循环
Python 的默认实现(称为 CPython)执行操作的速度非常慢。这部分是由于语言的动态解释性质:类型是灵活的,因此无法将操作序列编译为高效的机器代码,如 C 和 Fortran 等语言。
最近有各种解决这个弱点的尝试:众所周知的例子是 PyPy 项目,Python 的即时编译实现;Cython 项目,它将 Python 代码转换为可编译的 C 代码; 和 Numba 项目,它将 Python 代码片段转换为快速 LLVM 字节码。
每种方法都有其优点和缺点,但可以肯定的是,这三种方法都没有超过标准 CPython 引擎的范围和普及程度。Python 的相对迟缓通常体现在重复许多小操作的情况下 - 例如通过循环遍历数组来操作每个元素。
例如,假设我们有一个数组,我们想计算每个值的倒数。直截了当的方法可能如下所示:
代码语言:javascript复制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 背景的人来说,这种实现可能是相当自然的。但是如果我们对较大输入测量这个代码的执行时间,我们会发现这个操作非常慢,或许令人惊讶!
我们将使用 IPython 的%timeit
魔术指令(在“代码的性能度量和计时”中讨论)对此进行基准测试:
big_array = np.random.randint(1, 100, size=1000000)
%timeit compute_reciprocals(big_array)
# 1 loop, best of 3: 2.91 s per loop
计算这些数百万次操作并存储结果需要几秒钟!甚至当手机的处理速度以千兆 FLOPS 测量(即每秒数十亿次数值运算)时,这看起来几乎是非常缓慢的。
事实证明,这里的瓶颈不是操作本身,而是 CPython 必须在循环的每个循环中执行的类型检查和函数调度。每次计算倒数时,Python 首先检查对象的类型,并动态查找要用于该类型的正确函数。
如果我们使用编译代码,那么在代码执行之前就会知道这种类型规范,并且可以更有效地计算结果。
UFuncs 简介
对于许多类型的操作,NumPy 为这种静态类型的编译例程提供了方便的接口。 这称为向量化操作。实现方式为,简单地对数组执行操作,然后将该操作应用于每个元素。这种向量化方法旨在将循环推入 NumPy 背后的编译层,从而加快执行速度。
比较以下两个结果:
代码语言:javascript复制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 ]
'''
查看我们的大型数组的执行时间,我们发现它比 Python 循环快了几个数量级:
代码语言:javascript复制%timeit (1.0 / big_array)
# 100 loops, best of 3: 4.6 ms per loop
NumPy 中的向量化操作是通过ufunc
实现的,其主要目的是,对 NumPy 数组中的值快速执行重复操作。ufunc
非常灵活 - 在我们看到标量和数组之间的操作之前,我们也可以在两个数组之间操作:
np.arange(5) / np.arange(1, 6)
# array([ 0. , 0.5 , 0.66666667, 0.75 , 0.8 ])
ufunc
操作不仅限于一维数组 - 它们也可以作用于多维数组:
x = np.arange(9).reshape((3, 3))
2 ** x
'''
array([[ 1, 2, 4],
[ 8, 16, 32],
[ 64, 128, 256]])
'''
使用ufunc
向量化的计算,几乎总是比使用 Python 循环实现的对应方案更有效,特别是当数组的大小增加时。
每次在 Python 脚本中看到这样的循环时,都应该考虑是否可以用向量化表达式替换它。
探索 NumPy ufunc
ufunc
有两种形式:一元ufunc
,它在单个输入上运行,二元ufunc
,在两个输入上运行。我们将在这里看到这两种函数的例子。
数组算数
NumPy 的ufunc
使用起来非常自然,因为它们使用了 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]
'''
还有一个用于取负的一元ufunc
,一个用于求幂的**
运算符,以及一个用于取模的%
运算符:
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]
'''
此外,这些可以按照你的意愿串联在一起,并且遵守标准运算顺序:
代码语言:javascript复制-(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 | 描述 |
---|---|---|
| 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) |
另外还有布尔/位运算符; 我们将在“比较,掩码和布尔逻辑”中探索这些内容。
绝对值
就像 NumPy 理解 Python 的内置算术运算符一样,它也理解 Python 内置的绝对值函数:
代码语言:javascript复制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.])
三角函数
NumPy 提供了大量有用的ufunc
,对数据科学家来说最有用的是三角函数。我们首先定义一个角度数组:
theta = np.linspace(0, np.pi, 3)
现在我们可以在这些值上计算一些三角函数:
代码语言:javascript复制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.00000000e 00 1.00000000e 00 1.22464680e-16]
cos(theta) = [ 1.00000000e 00 6.12323400e-17 -1.00000000e 00]
tan(theta) = [ 0.00000000e 00 1.63312394e 16 -1.22464680e-16]
'''
这些值是在机器精度内计算的,这就是为什么,应该为零的值并不总是精确为零。反三角函数也可用:
代码语言:javascript复制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]
'''
指数和对数
NumPy ufunc
中另一种常见的操作类型是指数:
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. ]
'''
还有一些专用版本可用于为非常小的输入保持精度:
代码语言:javascript复制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
更精确的值。
专用ufunc
NumPy 还有更多的ufunc
可用,包括双曲线三角函数,按位算术,比较运算符,从弧度到度数的转换,舍入和余数等等。浏览 NumPy 文档可以发现许多有趣的函数。
ufunc
的另一个优秀来源是子模块scipy.special
,更专业和更隐蔽。如果你想对你的数据计算一些不常见的数学函数,它们很可能在scipy.special
中实现。
这些函数太多了,难以列出,但是下面的代码段显示了一些东西,可能出现在统计信息上下文中:
代码语言:javascript复制from scipy import special
# Gamma 函数(广义阶乘)和相关函数
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.00000000e 00 2.40000000e 01 3.62880000e 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]
'''
NumPy 和scipy.special
中有更多可用的ufunc
。由于这些软件包的文档可在线获取,因此搜索gamma function python
通常会找到相关信息。
高级ufunc
特性
许多 NumPy 用户在没有学习完整特性的情况下使用ufunc
。我们将在这里概述ufunc
的一些专用特性。
指定输出
对于大型计算,指定存储计算结果的数组,有时很有用。它不会创建临时数组,可以用于将计算结果直接写入你希望的内存位置。对于所有ufunc
,可以使用函数的out
参数来完成:
x = np.arange(5)
y = np.empty(5)
np.multiply(x, 10, out=y)
print(y)
# [ 0. 10. 20. 30. 40.]
这甚至可以用于数组视图。例如,我们可以将计算结果写入指定数组的每个其他元素:
代码语言:javascript复制y = np.zeros(10)
np.power(2, x, out=y[::2])
print(y)
# [ 1. 0. 2. 0. 4. 0. 8. 0. 16. 0.]
如果我们改为编写y [:: 2] = 2 ** x
,这将创建一个临时数组来保存2 ** x
的结果,然后将这些值复制到y
数组中。
对于如此小的计算而言,这并没有多大区别,但对于非常大的数组,通过小心使用out
参数可以节省大量内存。
聚合
对于二元ufunc
,有一些有趣的聚合可以从对象直接计算。例如,如果我们想要使用特定操作简化数组,我们可以使用任何ufunc
的reduce
方法。
reduce
会重复将给定操作应用于数组元素,直到只剩下一个结果。例如,在add ufunc
上调用reduce
会返回数组中所有元素的总和:
x = np.arange(1, 6)
np.add.reduce(x)
# 15
类似地,在multiply
ufunc上调用reduce
会产生所有数组元素的乘积:
np.multiply.reduce(x)
# 120
如果我们想存储计算的所有中间结果,我们可以使用accumulate
:
np.add.accumulate(x)
# array([ 1, 3, 6, 10, 15])
np.multiply.accumulate(x)
# array([ 1, 2, 6, 24, 120])
请注意,对于这些特殊情况,有专门的 NumPy 函数来计算结果(np.sum
,np.prod
,np.cumsum
,np.cumprod
), 我们将在“聚合:最小、最大和之间的任何东西”中探索。
外积
最后,任何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
方法也非常有用,我们将在“花式索引”中探索。
ufunc
的另一个非常有用的功能是,能够在不同大小和形状的数组之间操作,称为“广播”。这个主题非常重要,我们将为它编写一整节(参见“数组计算:广播”)。
ufunc
:了解更多
通用函数的更多信息(包括可用函数的完整列表)可在 NumPy 和 SciPy 文档站点上找到。
回想一下,你也可以通过导入软件包,并使用 IPython 的 TAB 补全和帮助(?
)功能,直接从 IPython 中访问信息,如“IPython 中的帮助和文档”中所述。