NumPy 数组计算：通用方法

NumPy 数组计算：通用方法

迄今为止，我们已探讨了NumPy的基础知识。接下来的章节中，我们将深入解析NumPy在Python数据科学领域占据重要地位的根本原因：其通过提供简洁灵活的接口，可实现对数据数组的高效计算。

NumPy数组的计算可能极其高效，也可能异常缓慢。
实现高效计算的关键在于采用向量化操作，这些操作通常通过NumPy的通用函数（universal functions，简称ufuncs）来实现。
本章将深入阐述NumPy通用函数的必要性——它们能显著提升对数组元素进行重复计算的效率，并系统介绍NumPy包中最常用且实用的算术通用函数。

一、循环的低效性

Python默认实现（即CPython）在执行某些操作时效率较低。这主要源于其动态解释型语言的特性：类型系统高度灵活，导致操作序列无法像C、Fortran等语言那样被编译为高效的机器码。近年来业界提出了多种改进方案，典型代表包括：即时编译型Python实现PyPy项目；可将Python代码转换为可编译C代码的Cython项目；以及能将Python代码片段转换为快速LLVM字节码的Numba项目。尽管这些方案各有优劣，但可以确定的是，目前尚未有任何一种方法能完全取代标准CPython引擎的生态影响力和应用普及度。

Python的相对性能迟滞通常集中体现在需要重复执行大量小型操作的场景中，具体而言，例如遍历数组并对每个元素进行操作的情况。假设我们有一个数值数组，需要计算其中每个元素的倒数，一种直观的实现方式可能如下所示：

import numpy as np
rng = np.random.default_rng(seed=1701)def compute_reciprocals(values):output = np.empty(len(values))for i in range(len(values)):output[i] = 1.0 / values[i]return outputvalues = rng.integers(1, 10, size=5)
compute_reciprocals(values)

对于有C或Java背景的开发者而言，这种编程方式可能显得相当自然。但当我们测量该代码处理大规模输入时的执行时间，就会发现这种操作极其低效——其缓慢程度可能令人咋舌！我们将使用IPython的%timeit魔法命令进行基准测试。

big_array = rng.integers(1, 100, size=1000000)
%timeit compute_reciprocals(big_array)

完成百万次运算并存储结果竟需数秒（原文，实测本机约1秒）！在移动设备都已具备千兆次浮点运算（gigaflops，即每秒数十亿次数值运算）能力的今天，这种效率低下几乎令人难以置信。事实上，瓶颈并非运算本身，而在于CPython在每次循环中必须执行的类型检查与函数派发（function dispatches）机制——每当计算倒数时，Python解释器首先需要动态解析操作数类型，并查找适用于该类型的正确函数。若采用编译型语言实现，类型信息可在代码执行前确定，从而显著提升计算效率。

二、通用函数（Ufuncs）介绍

NumPy为各类运算提供了便捷的接口，可直接对接这种静态类型的编译层例程。这种机制被称为向量化操作。对于诸如逐元素除法这类简单运算，向量化的实现方式就如同直接在数组对象上使用Python算术运算符般直观。这种设计理念的本质在于将循环操作下沉至NumPy底层的编译层，从而实现执行效率的飞跃式提升。

请比较下列两种操作的执行效能差异：

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

[0.11111111 0.25       1.         0.33333333 0.125     ]
[0.11111111 0.25       1.         0.33333333 0.125     ]

观察大规模数组的执行时间，可以发现其完成速度较Python循环提升达数个数量级。

%timeit (1.0 / big_array)

2.68 ms ± 18.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

NumPy的向量化操作通过通用函数（ufuncs）实现，其核心设计目标是对数组元素进行高效重复运算。通用函数具有高度扩展性——此前的示例展示了标量与数组间的运算，而更可实现双数组间的元素级操作：

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

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

通用函数的运算维度并不仅限于一维数组，其更可实现对多维数组的跨维度操作：

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

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

通过通用函数(ufuncs)实现的向量化计算，始终具备较Python循环方案更显著的效率优势。当数组规模持续增大时，这种效率优势呈指数级扩大。在NumPy编程实践中，每当遇到循环结构时，开发者应当优先考量能否将其重构为向量化表达式——这已成为高性能科学计算的最佳实践准则。

三、NumPy通用函数深度解析

通用函数(ufuncs)主要分为两类：一元通用函数，作用于单一输入；二元通用函数，处理双输入操作。本章将通过具体案例解析这两类函数的应用场景及实现机制。

1. 数组运算（算术）

NumPy通用函数的使用体验高度符合直觉，这得益于其完美复用了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)   # floor division 往下取整

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]

此外还包括三类特殊运算符：用于数值取反的一元通用函数（unary 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内置特定通用函数（ufuncs）的语法糖。以加法运算符+为例，其底层实现实际上是对应着add通用函数的具体调用：

np.add(x, 2)

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

下表列明了NumPy实现的算术运算符及其对应机制：

运算符	等效通用函数	功能描述
+	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还实现了布尔逻辑与位运算两类运算符体系。

2. 绝对值

正如NumPy原生支持Python内置的算术运算符，其同样兼容处理Python绝对值函数（abs()）：

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

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

该通用函数亦可处理复数型数据，此时其将返回复数的模（即复平面上的向量长度）：

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

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

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]

计算结果在机器精度范围内，这解释了理论上应为零的数值可能不会精确为零的现象。NumPy同样提供了反三角函数集：

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三角函数相关通用函数小结

NumPy提供完整的三角函数集，支持以弧度为单位的向量化计算：

基础三角运算

函数	描述	数学表达式
np.sin	正弦函数	sinθ
np.cos	余弦函数	cosθ
np.tan	正切函数	tanθ = sinθ/cosθ

反三角函数

函数	描述	值域
np.arcsin	反正弦函数	[-π/2, π/2]
np.arccos	反余弦函数	[0, π]
np.arctan	反正切函数	(-π/2, π/2)

角度转换

函数	转换方向
np.deg2rad	角度 → 弧度
np.rad2deg	弧度 → 角度

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

作为指数函数的逆运算，对数函数同样在NumPy通用函数体系中得以实现。基础函数np.log对应自然对数运算；若需计算以2为底或以10为底的对数，NumPy亦提供专门实现：

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

针对极小输入值，NumPy还提供若干专用函数变体，这些函数能有效保持计算精度：

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函数，这些专用函数可输出更高精度的计算结果（因避免舍入误差）。

技术原理说明：

• 常规运算exp(x)-1在x趋近0时，由于浮点数的有限精度（IEEE 754双精度约16位有效数字），计算结果会被舍入为0
• expm1采用泰勒展开近似计算：当|x| < 1e-5时，expm1(x) ≈ x + x²/2 + x³/6（保留更多有效位数）
• log1p通过重构运算路径避免1+x的精度损失，直接计算ln(1+x)的泰勒展开式

指数函数和对数函数内容小结

指数函数

函数	数学表达式	功能描述
np.exp	$e^x$	自然指数函数
np.exp2	$2^x$	基数2指数函数
np.expm1	$e^x-1$	数值稳定的指数减一运算

对数函数

函数	数学表达式	功能描述
np.log	$ln(x)$	自然对数函数
np.log10	$lo{g}_{10}(x)$	常用对数（基10）
np.log2	$lo{g}_2(x)$	二进制对数（基2）
np.log1p	$ln(1+x)$	数值稳定的对数运算

5. 专用通用函数

除前述基础功能外，NumPy还提供了更为丰富的通用函数体系，涵盖双曲三角函数、位级运算、比较运算、弧度与角度的相互转换、舍入与求余运算等众多领域。深入研读NumPy官方文档可发掘大量高阶数学功能。

对于更为专业化的通用函数需求，推荐使用scipy.special子模块。该模块几乎囊括了所有已知的数学函数实现——若需对数据执行某些冷僻的数学变换，极可能已在scipy.special中实现。由于函数数量过于庞大，此处仅例示统计学中常用的几种函数：

from scipy import special

# Gamma functions (generalized factorials) and related functions
# Gamma 函数 （广义阶乘）和相关函数
x = [1, 5, 10]
print("gamma(x)     =", special.gamma(x))   # gamma(x) = (x-1)!
print("ln|gamma(x)| =", special.gammaln(x)) # ln|gamma(x)|
print("beta(x, 2)   =", special.beta(x, 2)) # beta(x, y) = gamma(x) * gamma(y) / gamma(x + y)

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]

# 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模块中集成的通用函数（ufuncs）体系极为庞大，其规模远超基础教程所能涵盖。得益于完善的在线文档支持，开发者通过"gamma函数 Python"等关键词进行网络搜索，即可精准定位相关函数的官方技术文档。

四、通用函数高阶特性

多数NumPy使用者虽常运用通用函数(ufuncs)，却未必深究其完整功能特性。本文将重点解析通用函数的若干专用特性。

1. 指定输出目标

对于大规模数值计算，通过指定结果存储数组可有效优化内存使用。所有通用函数(ufuncs)均支持通过out参数实现输出目标预设：

x = np.arange(5)
y = np.empty(5)# 预分配输出数组
np.multiply(x, 10, out=y)
print(y)  # 输出: [ 0. 10. 20. 30. 40.]# 支持视图操作
y_center = y[1:-1]
np.power(y_center, 2, out=y_center) 
print(y)  # 输出: [ 0. 100. 400. 900. 40.]

[ 0. 10. 20. 30. 40.]
[  0. 100. 400. 900.  40.]

该特性甚至可与数组视图结合使用。例如，我们可以将计算结果写入指定数组的间隔元素位置：

y = np.zeros(10)
np.power(2, x, out=y[::2])  # 每隔一个元素赋值
print(y)  # 输出: [ 1. 0. 2. 0. 4. 0. 8. 0. 16. 0.]

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

若采用常规写法y[::2] = 2 ** x，系统将经历两个内存密集型操作：一是临时创建数组用于存储2 ** x的计算结果，二是接下来的数据复制操作，即将计算结果所在的临时数组复制到数组y中。虽然在小规模计算中np.power和**两中求幂方式差异甚微，但在处理超大规模数组时，通过合理运用out参数的精细化内存管理，可节省更多内存。

2. 聚合操作

对于二元通用函数(ufuncs)，可直接通过其对象实现聚合计算。以归约运算为例，任何ufunc均可调用reduce方法对数组执行持续聚合操作——该方法将持续应用特定运算直至生成单一结果。

例如，在加法通用函数add上调用reduce，可获得数组所有元素的累加和：

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

类似地，在乘法通用函数multiply上调用reduce方法，将生成数组所有元素连乘的积：

np.multiply.reduce(x)   # 120

若需保留计算过程中各阶段的中间结果，可转而调用accumulate方法。该方法将生成一个与输入数组等长的输出数组，其中每个元素代表截止当前索引位置的累积运算结果：

np.add.accumulate(x)    # [ 1  3  6 10 15]

np.multiply.accumulate(x)   # [  1   2   6  24 120]

需要特别指出的是，针对前文所述的归约与累积运算场景，NumPy已提供经过深度优化的专用函数实现（如np.sum、np.prod、np.cumsum、np.cumprod）。

注：实际运行效率未必时最佳（需要实际代码验证），但这些专用函数确实便于使用。

# import numpy as np
arr = np.array([1, 2, 3, 4])# 性能对比测试
%timeit np.add.reduce(arr)   
%timeit np.sum(arr)          %timeit np.multiply.accumulate(arr)  
%timeit np.cumprod(arr)      

512 ns ± 3.02 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
1.05 μs ± 11.5 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
287 ns ± 3.49 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
771 ns ± 5.4 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

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)在不同形状与尺寸数组间的运算能力——这一系列智能维度扩展操作被称作广播（broadcasting）。

# import numpy as np# 典型广播案例
A = np.array([[1, 2, 3],          # 形状 (2,3)[4, 5, 6]])
B = np.array([10, 20, 30])       # 形状 (3,)# 广播机制自动扩展B到(2,3)后执行运算
print(A + B)
"""
输出:
[[11 22 33][14 25 36]]
"""

关于NumPy通用函数的高级特习惯，可总结如下，供参考：

五、通用函数：深度探索指南

有关通用函数（包含完整函数列表）的更多技术细节，请参阅以下权威文档源：

• ​​NumPy官方文档​​ - 提供API参考与最佳实践指南；
• ​​SciPy技术文档​​ - 涵盖高级科学计算函数库。

需特别强调的是，开发者可在IPython环境中通过以下工作流直接获取技术信息：

• 1. ​导入目标包​​：加载NumPy/SciPy等计算库；
• 2. ​​使用Tab自动补全​​：探索函数/属性列表；
• 3. ​​调用帮助系统​​：通过?操作符查询文档细节。

示例代码如下：

# 操作示例
In [1]: import numpy as np# Tab补全探索（输入np.后按Tab）
In [2]: np.[Tab]
np.abs    np.add     np.array   ...（显示所有以np.开头的函数）# 帮助查询（输入函数名加?）
In [3]: np.abs?
Signature: np.abs(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])
Docstring:
Calculate the absolute value element-wise.
...（显示完整文档字符串）# 即时性能测试
In [4]: %timeit np.abs(np.random.randn(1000))
14.1 µs ± 437 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)