9.5 NumPy 数组上的计算: 通用函数
本节是《Python 数据科学手册》(Python Data Science Handbook)的摘录.
译者: 飞龙 https://github.com/wizardforcel
协议: CC BY-NC-SA 4.0 http://creativecommons.org/licenses/by-nc-sa/4.0/
到目前为止, 我们一直在讨论 NumPy 的一些基本要点; 在接下来的几节中, 我们将深入探讨 NumPy 在 Python 数据科学领域如此重要的原因.
也就是说, 它为数据数组的最优计算, 提供了一个简单而灵活的接口.
NumPy 数组的计算速度非常快, 也可能非常慢. 使其快速的关键是使用向量化操作, 通常通过 NumPy 的通用函数 (ufunc) 实现.
本节激发了 NumPy 的 ufunc 的需求, 这些 ufunc 可用于更有效地对数组元素进行重复计算. 然后介绍了 NumPy 包中可用的, 许多最常用和最有用的算术 ufunc.
慢速的循环
Python 的默认实现 (称为 CPython) 执行操作的速度非常慢.
这部分是由于语言的动态解释性质: 类型是灵活的, 因此无法将操作序列编译为高效的机器代码, 如 C 和 Fortran 等语言.
最近有各种解决这个弱点的尝试: 众所周知的例子是 PyPy http://pypy.org/ 项目, Python 的即时编译实现; Cython http://cython.org 项目, 它将 Python 代码转换为可编译的 C 代码; 和 Numba http://numba.pydata.org/ 项目, 它将 Python 代码片段转换为快速 LLVM 字节码.
每种方法都有其优点和缺点, 但可以肯定的是, 这三种方法都没有超过标准 CPython 引擎的范围和普及程度.
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 背景的人来说, 这种实现可能是相当自然的.
但是如果我们对较大输入测量这个代码的执行时间, 我们会发现这个操作非常慢, 或许令人惊讶!
我们将使用 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 背后的编译层, 从而加快执行速度.
比较以下两个结果:
- 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 循环快了几个数量级:
- %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]
- '''
此外, 这些可以按照你的意愿串联在一起, 并且遵守标准运算顺序:
- -(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 内置的绝对值函数:
- 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)
现在我们可以在这些值上计算一些三角函数:
- 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]
- '''
这些值是在机器精度内计算的, 这就是为什么, 应该为零的值并不总是精确为零. 反三角函数也可用:
- 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. ]
- '''
还有一些专用版本可用于为非常小的输入保持精度:
- 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 中实现.
这些函数太多了, 难以列出, 但是下面的代码段显示了一些东西, 可能出现在统计信息上下文中:
- 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.]
这甚至可以用于数组视图. 例如, 我们可以将计算结果写入指定数组的每个其他元素:
- 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 http://www.numpy.org 和 SciPy http://www.scipy.org 文档站点上找到.
回想一下, 你也可以通过导入软件包, 并使用 IPython 的 TAB 补全和帮助 (?) 功能, 直接从 IPython 中访问信息, 如 "IPython 中的帮助和文档" 中所述.
来源: https://juejin.im/entry/5c2ef2b8518825261e1eea68