< Numpy数组基础 | 目录 | 聚合:Min, Max, 以及其他 >
直到目前为止,我们已经讨论了一些NumPy的基本构件;在下面几个小节中,我们会深入讨论NumPy能在Python数据科学中占据重要地位的原因。简而言之,NumPy提供了简单和灵活的接口来对数组数据计算进行优化。
对NumPy的数组进行计算相较其他普通的实现方式而言是非常快的。快的原因关键在于使用了向量化的操作,因为它们都是通过NumPy的通用函数(ufuncs)实现的。希望通过本节的介绍,能让读者习惯使用ufuncs,它们能使在数组元素上的重复计算更加快速和高效。本节还会介绍许多NumPy中最常用的ufuncs数学计算方法。
循环,慢的实现
Python的默认实现(被称为CPython)对于一些操作执行效率很低。这部分归咎于语言本身的动态和解释执行特性:因为类型是动态的,因此不到执行时,无法预知变量的类型,因此不能像C或者Fortran那样预先将代码编译成机器代码来执行。近年来,也出现了很多尝试来弥补这个缺陷:其中比较流行和著名的包括PyPy,Python的JIT编译实现;Cython,可以将Python代码转换为可编译的C代码;和Numba,可以将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)
1.63 s ± 3.21 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
这个操作对于百万级的数据集耗时需要几秒。当现在手机的每秒浮点数运算次数都已经已经达到10亿级别,这实在是不可思议的慢了。通过分析发现瓶颈并不是代码本身,而是每次循环时CPython必须执行的类型检查和函数匹配。每次计算倒数时,Python首先需要检查对象的类型,然后寻找一个最合适的函数对这种类型进行计算。如果我们使用编译型的语言实现上面的代码,每次计算的时候,类型和应该执行的函数都已经确定,因此执行的时间肯定短很多。
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)
1.48 ms ± 53.2 µs per loop (mean ± std. dev. of 7 runs, 1,000 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]])
通过ufuncs向量化计算基本上都会比使用Python循环实现的相同方法要更加高效,特别是数组的长度增长的情况下。任何情况下,如果你看到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 ) |
除此之外还有布尔和二进制位操作;我们会在比较,遮盖和布尔逻辑中介绍它们。
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.])
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]
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
计算得到更加精确的结果。
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"就可以找到相关的信息。
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
参数能节省很多内存空间。
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
),我们会在聚合:Min, Max, 以及其他中详细讨论。
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]])
更多有关ufuncs的信息(包括完整的函数列表)可以在NumPy 和 SciPy的在线文档获得。
不要忘记了我们可以使用IPython的帮助工具?
来获取任何相关的帮助信息,正如我们在IPython的帮助和文档中介绍过的那样。
< Numpy数组基础 | 目录 | 聚合:Min, Max, 以及其他 >