Numpy介绍(3)数组通用函数

NumPy 数组的计算可以非常快,也可以非常慢。使其快速的关键是使用矢量化操作,通常通过 NumPy 的通用函数 (ufuncs) 实现。
本节激发了对 NumPy 的 ufunc 的需求,它可用于更有效地对数组元素进行重复计算。然后介绍了 NumPy 包中提供的许多最常见和有用的算术ufunc。

循环的缓慢

Python 的默认实现(称为 CPython)执行某些操作的速度非常慢。这在一定程度上是由于语言的动态、解释性质:类型是灵活的,因此操作序列不能像 C 和 Fortran 等语言那样被编译成高效的机器代码。最近有各种尝试来解决这个弱点:众所周知的例子是 PyPy 项目,它是 Python 的即时编译实现; Cython 项目,它将 Python 代码转换为可编译的 C 代码;和 Numba 项目,它将 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 magic函数对其进行基准测试:

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

1 loop, best of 3: 2.91 s per loop

计算这百万次操作并存储结果需要几秒钟!当手机的处理速度以 Giga-FLOPS(即每秒数十亿次数值运算)衡量时,这似乎慢得几乎荒谬。事实证明,这里的瓶颈不是操作本身,而是 CPython 必须在循环的每个周期执行的类型检查和函数调度。每次计算倒数时,Python 首先检查对象的类型并动态查找用于该类型的正确函数。如果我们在编译代码中工作,那么在代码执行之前就会知道这种类型规范,并且可以更有效地计算结果。

介绍UFuncs

对于许多类型的操作,NumPy 为这种静态类型的编译例程提供了一个方便的接口。这称为矢量化操作。这可以通过简单地对数组执行操作来完成,然后将其应用于每个元素。这种矢量化方法旨在将循环推入作为 NumPy 基础的编译层,从而加快执行速度。
比较一下下面2个的结果:

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 中的矢量化操作是通过 ufuncs 实现的,其主要目的是快速对 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的UFuncs

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])

下表列出了 NumP 中实现的算术运算符:

操作符 等效的ufunc 描述
+ np.add 加 (e.g., 1 + 1 = 2)
- np.subtract 减 (e.g., 3 - 2 = 1)
- np.negative 一元取负 (e.g., -2)
/ np.divide 除 (e.g., 3 / 2 = 1.5)
// np.floor_divide 地板除 (e.g., 3 // 2 = 1)
** np.power 取幂 (e.g., 2 ** 3 = 8)
% np.mod 模数/余数 (e.g., 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 时更精确。

专门的ufuncs

NumPy 有更多可用的 ufunc,包括双曲三角函数、按位算术、比较运算符、从弧度到度的转换、舍入和余数等等。浏览 NumPy 文档会发现很多有趣的功能。
更专业和晦涩的 ufunc 的另一个极好的来源是子模块 scipy.special。如果你想对你的数据计算一些晦涩的数学函数,它很可能是在 scipy.special 中实现的。函数太多无法一一列举,但以下代码片段显示了统计上下文中可能出现的几个函数:

from scipy import special
# Gamma functions (generalized factorials) and related functions
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]
# Error function (integral of Gaussian)
# its complement, and its inverse
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”这样的网络搜索一般都能找到相关资料。

进阶的Ufuncs属性

我们将在这里概述 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

同样,对乘法 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 的另一个非常有用的特性是能够在不同大小和形状的数组之间进行操作,一组操作称为广播(broadcasting)。这个主题很重要,我们将用一整节来讨论它。

Ufuncs:更多探索

有关通用函数的更多信息(包括可用函数的完整列表)可以在 NumPy
SciPy 文档网站上找到。
回想一下,您还可以通过导入包并使用 IPython 的制表符完成和帮助 (?) 功能直接从 IPython 中访问信息,如Help and Documentation in IPython
中所述。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,504评论 6 496
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,434评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 160,089评论 0 349
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,378评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,472评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,506评论 1 292
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,519评论 3 413
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,292评论 0 270
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,738评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,022评论 2 329
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,194评论 1 342
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,873评论 5 338
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,536评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,162评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,413评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,075评论 2 365
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,080评论 2 352

推荐阅读更多精彩内容