Python 基础体系 · 第 94/112 篇。示例统一以 Python 3.14 为语言基线；第三方库使用与其兼容的现代稳定版本，版本敏感行为会单独说明。

NumPy 完整基础：ndarray、形状、广播、索引、向量化和内存布局

NumPy 是 Python 科学计算生态的底层数组库，核心对象是多维数组 ndarray，并围绕它提供数据类型、形状变换、索引、广播、通用函数、统计、线性代数、随机数和输入输出等能力。官方文档将 ndarray、广播、索引、数据类型、视图与副本、ufunc 等列为基础主题。(numpy.org)

本文以 Python 3.14 为代码运行范围，重点解释 NumPy 的数据模型和执行机制。NumPy 本身由 C、C++、Fortran 等底层实现支撑，但使用者首先需要理解：一个数组不是“嵌套 Python 列表的快捷版本”，而是一块由数据缓冲区、数据类型、形状和步长共同描述的内存视图。

一、安装与最小运行环境

建议在虚拟环境中安装 NumPy，避免不同项目共享依赖导致版本和二进制扩展冲突。官方安装文档同时提供 venv、pip、conda、uv 等方式，并建议安装后打印版本进行验证。(numpy.org)

python3.14 -m venv .venv
source .venv/bin/activate        # macOS/Linux
# .venv\Scripts\activate         # Windows

python -m pip install numpy
python -c "import numpy as np; print(np.__version__)"

在 Windows 上，python3.14 可能对应为 py -3.14：

py -3.14 -m venv .venv
.venv\Scripts\activate
python -m pip install numpy

代码中通常使用别名：

import numpy as np

np 只是约定俗成的别名，不是 NumPy 的特殊语法。真正重要的是保持数组类型来自同一个 NumPy 环境，并通过 np.__version__、np.show_config() 等方式排查运行环境。

二、ndarray 的本质：数据、类型、形状和步长

2.1 ndarray 不是普通列表

Python 列表保存的是对象引用：

values = [1, 2, 3]

列表中的每个元素可以是不同类型，列表本身只负责管理引用。NumPy 的 ndarray 通常保存一段具有统一数据类型的连续或规则间隔内存：

a = np.array([1, 2, 3], dtype=np.int64)

可以把一个数组抽象成：

A=(data,dtype,shape,strides)A = (\text{data}, \text{dtype}, \text{shape}, \text{strides})

其中：

  • data：底层数据缓冲区；
  • dtype：每个元素的解释方式和字节大小；
  • shape：每个轴上的元素数量；
  • strides：沿每个轴移动一个元素时需要跳过多少字节。

例如：

a = np.array([[10, 20, 30],
              [40, 50, 60]], dtype=np.int64)

print(a.ndim)      # 2
print(a.shape)     # (2, 3)
print(a.size)      # 6
print(a.dtype)     # int64
print(a.itemsize)  # 8
print(a.strides)   # 通常为 (24, 8)

shape == (2, 3) 表示数组有两个轴，第一轴长度为 2，第二轴长度为 3。由于 int64 每个元素占 8 字节，一行有 3 个元素，因此沿第 0 轴移动一格通常需要跳过：

3×8=24 字节3 \times 8 = 24\text{ 字节}

沿第 1 轴移动一格只需跳过 8 字节。

strides 是理解转置、切片、连续性和性能的关键，而不是一个只用于调试的附属属性。

2.2 ndim、shape、size 和 itemsize

x = np.zeros((2, 3, 4), dtype=np.float32)

print(x.ndim)      # 3
print(x.shape)     # (2, 3, 4)
print(x.size)      # 24
print(x.itemsize)  # 4
print(x.nbytes)    # 96

这些属性分别表示：

属性 含义
ndim 轴的数量，也就是维度数
shape 每个轴的长度
size 所有元素总数，等于各轴长度的乘积
itemsize 单个元素占用的字节数
nbytes 数组元素数据占用的总字节数，通常等于 size * itemsize

对于 shape = (2, 3, 4)：

size=2×3×4=24\text{size} = 2 \times 3 \times 4 = 24

注意，nbytes 只描述元素缓冲区大小，不一定包括 Python 对象本身、对象型数组中被引用对象的实际内存，以及数组视图的额外元数据。

2.3 轴不是“方向”，而是一个索引维度

考虑二维数组：

a = np.array([
    [1, 2, 3],
    [4, 5, 6],
])

它的形状是 (2, 3)：

  • 第 0 轴长度为 2，对应两行；
  • 第 1 轴长度为 3，对应三列。
print(a.sum(axis=0))  # [5 7 9]
print(a.sum(axis=1))  # [ 6 15]

axis=0 表示沿第 0 轴聚合。第 0 轴上的元素被压缩，因此保留第 1 轴：

[123456]⟶[1+42+53+6]\begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix} \longrightarrow \begin{bmatrix} 1+4 & 2+5 & 3+6 \end{bmatrix}

axis=1 则压缩列这个轴，得到每一行的和：

[1+2+34+5+6]\begin{bmatrix} 1+2+3\\ 4+5+6 \end{bmatrix}

使用 keepdims=True 可以保留被聚合轴的长度 1：

print(a.sum(axis=1, keepdims=True))
# [[ 6]
#  [15]]

这对于后续广播很有用，因为结果形状从 (2,) 变成了 (2, 1)。

三、创建数组与数据类型

3.1 从 Python 数据创建

a = np.array([1, 2, 3])
b = np.array([[1, 2], [3, 4]])

嵌套序列必须能够形成规则形状：

np.array([[1, 2], [3, 4]])       # shape 为 (2, 2)

如果嵌套层级或长度不一致，现代 NumPy 通常会拒绝直接创建规则数值数组：

# 不应依赖这种写法：
# np.array([[1, 2], [3]])

不规则嵌套数据本质上不是常规的多维数值数组。即使显式指定 dtype=object 使其能够保存，也会失去大量 NumPy 数值计算优势：

ragged = np.array([[1, 2], [3]], dtype=object)

对象型数组保存的是 Python 对象引用，运算往往回到 Python 层，不能等同于固定宽度数值数组。

3.2 全零、全一、单位矩阵和范围数组

np.zeros((2, 3))
np.ones((2, 3), dtype=np.int32)
np.empty((2, 3))
np.eye(3)
np.arange(0, 10, 2)
np.linspace(0, 1, 5)

区别如下：

  • zeros：元素初始化为零；
  • ones：元素初始化为一；
  • empty：只分配内存，不初始化元素；
  • eye：创建单位矩阵；
  • arange：按步长生成，停止值通常不包含在内；
  • linspace：按元素数量在区间内均匀采样。

empty 的值不可预测：

x = np.empty(5)
print(x)

它可能显示之前该内存区域遗留的任意位模式。只有在后续会完整覆盖所有位置时，empty 才有意义：

x = np.empty(5)
x[:] = 0

对于浮点步长，通常优先考虑 linspace：

np.arange(0, 1, 0.1)
np.linspace(0, 1, 11)

arange 根据步长不断累加，浮点舍入可能使元素数量和边界表现不符合直觉；linspace 直接根据目标数量计算采样点。

3.3 dtype：数组如何解释位模式

a = np.array([1, 2, 3], dtype=np.int32)
b = np.array([1, 2, 3], dtype=np.float64)
c = np.array([True, False], dtype=np.bool_)

常见数值类型包括：

  • bool：布尔值；
  • int8、int16、int32、int64：有符号整数；
  • uint8、uint16、uint32、uint64：无符号整数；
  • float16、float32、float64：浮点数；
  • complex64、complex128：复数；
  • 字符串、字节串和结构化数据类型。

类型不是装饰信息，它决定：

  1. 每个元素占用多少字节；
  2. 可表示的数值范围；
  3. 运算时的类型提升规则；
  4. 与外部 C、Fortran 或文件格式交换时如何解释内存。

整数溢出

固定宽度整数不会自动变成任意精度整数：

x = np.array([127], dtype=np.int8)
print(x + 1)

在特定运算路径中可能发生溢出或类型提升。生产代码不能假设 NumPy 整数像 Python int 一样无限扩展。

浮点精度

x = np.array([0.1, 0.2], dtype=np.float64)
print(x.sum() == 0.3)  # 不应依赖为 True

二进制浮点数通常无法精确表示十进制小数。比较浮点结果应根据问题尺度使用容差：

np.isclose(x.sum(), 0.3)
np.allclose(x, expected, rtol=1e-5, atol=1e-8)

atol 是绝对误差容忍度，rtol 是相对误差容忍度。数值接近零时，绝对容差尤其重要。

3.4 类型转换和类型提升

a = np.array([1, 2, 3], dtype=np.int32)
b = a.astype(np.float64)

print(a.dtype)  # int32
print(b.dtype)  # float64

astype 默认创建新数组，而不是原地修改原数组。

不同类型运算时，NumPy 会根据类型提升规则选择结果类型：

a = np.array([1, 2], dtype=np.int32)
b = np.array([0.5, 1.5], dtype=np.float64)

print((a + b).dtype)  # 通常为 float64

不要只看输入数组类型，还要检查中间结果：

result = a + b
print(result.dtype)

对于需要严格控制内存和精度的代码，可以显式指定类型，而不是依赖隐式转换。

四、形状：数组结构的代数约束

4.1 reshape 只改变解释方式

x = np.arange(6)
y = x.reshape(2, 3)

print(x)
# [0 1 2 3 4 5]

print(y)
# [[0 1 2]
#  [3 4 5]]

reshape 要求元素总数不变：

∏iold_shapei=∏jnew_shapej\prod_i \text{old\_shape}_i = \prod_j \text{new\_shape}_j

因此：

np.arange(6).reshape(2, 3)  # 合法
np.arange(6).reshape(3, 2)  # 合法
# np.arange(6).reshape(4, 2) # ValueError

可以使用一个 -1 让 NumPy 自动推导该维度：

x.reshape(2, -1)  # 推导为 (2, 3)
x.reshape(-1, 3)  # 推导为 (2, 3)

只能有一个 -1，因为多个未知维度无法由元素总数唯一确定。

4.2 reshape 可能返回视图，也可能返回副本

如果新的形状能够通过重新解释原始步长得到，reshape 可能返回视图；否则可能分配副本。不能仅凭函数名判断是否共享内存。

x = np.arange(6)
y = x.reshape(2, 3)

y[0, 0] = 100
print(x)
# [100   1   2   3   4   5]

可以通过 np.shares_memory 进行验证：

print(np.shares_memory(x, y))

shares_memory 比检查 base 更直接，因为复杂视图可能经过多层派生。

4.3 ravel、flatten 和 reshape(-1)

x = np.arange(6).reshape(2, 3)

a = x.ravel()
b = x.flatten()
c = x.reshape(-1)

典型区别：

  • ravel()：尽量返回视图，必要时才复制；
  • flatten()：总是返回副本；
  • reshape(-1)：尝试改变形状，可能返回视图，也可能复制。

如果修改扁平结果不应影响原数组，使用 flatten() 或显式 .copy()：

flat = x.ravel().copy()

4.4 转置与轴置换

二维数组的 .T 交换轴：

x = np.array([
    [1, 2, 3],
    [4, 5, 6],
])

print(x.T)
# [[1 4]
#  [2 5]
#  [3 6]]

对于高维数组，使用 transpose 或 moveaxis 更明确：

x = np.zeros((2, 3, 4))

y = x.transpose(2, 0, 1)
print(y.shape)  # (4, 2, 3)

z = np.moveaxis(x, 0, -1)
print(z.shape)  # (3, 4, 2)

.T 对高维数组会反转轴顺序：

x.T.shape  # (4, 3, 2)

在图像和深度学习数据中，HWC、CHW、NHWC、NCHW 的差异就是轴顺序差异。形状数值相同并不代表语义相同。

4.5 添加和删除长度为 1 的轴

x = np.array([1, 2, 3])
print(x.shape)  # (3,)

a = x[:, None]
b = x[None, :]

print(a.shape)  # (3, 1)
print(b.shape)  # (1, 3)

None 等价于 np.newaxis，用于在指定位置插入长度为 1 的轴。

print(np.expand_dims(x, axis=0).shape)  # (1, 3)
print(np.expand_dims(x, axis=1).shape)  # (3, 1)

删除长度为 1 的轴：

x = np.zeros((1, 3, 1))
print(np.squeeze(x).shape)  # (3,)

如果只允许删除某个指定轴，可以传入 axis：

np.squeeze(x, axis=0)

这样当该轴长度不是 1 时会报错，能够防止意外删除错误的维度。

五、广播：不复制数据的形状匹配规则

广播是 NumPy 允许不同形状数组参与逐元素运算的规则。它不是任意形状自动对齐，而是对形状从右向左逐维比较。

两个维度兼容，当且仅当：

  1. 两个维度相等；或
  2. 其中一个维度为 1；或
  3. 某个数组在该位置没有维度，相当于缺失维度。

比较形状时，从最后一个轴开始。

5.1 标量广播

x = np.array([1, 2, 3])
print(x + 10)
# [11 12 13]

标量可以看作形状为空的数组。它与任何形状兼容，因此 10 被概念上应用到每一个元素。

5.2 一维数组与二维数组

x = np.array([
    [1, 2, 3],
    [4, 5, 6],
])

bias = np.array([10, 20, 30])

print(x + bias)
# [[11 22 33]
#  [14 25 36]]

形状比较：

x:     (2, 3)
bias:    (3)

从右侧比较：

3 与 3：相等
2 与缺失：兼容

结果形状为 (2, 3)。

这个例子表示给每一列添加一个偏置。若想给每一行添加偏置，必须让偏置形状为 (2, 1)：

row_bias = np.array([100, 200])[:, None]

print(x + row_bias)
# [[101 102 103]
#  [204 205 206]]

如果直接写：

x + np.array([100, 200])

形状为 (2, 3) 和 (2,)，从右侧比较时：

3 与 2：不相等，且都不为 1

因此会抛出广播错误。

5.3 广播的完整推导

设：

a.shape == (8, 1, 6)
b.shape == (7, 6)

为比较方便，给 b 在左侧补一个长度为 1 的轴：

a: (8, 1, 6)
b: (1, 7, 6)

逐轴比较：

8 与 1：兼容，结果为 8
1 与 7：兼容，结果为 7
6 与 6：相等，结果为 6

所以结果形状为：

(8, 7, 6)

但下面的形状不能广播：

a.shape == (8, 2, 6)
b.shape == (7, 6)

补齐后：

a: (8, 2, 6)
b: (1, 7, 6)

第二轴比较 2 与 7，既不相等，也没有一个是 1，所以失败。

5.4 广播不等于复制

从概念上，bias.shape == (3,) 被扩展为 (2, 3)，但 NumPy 通常不需要真的创建两行副本。它可以通过步长为 0 的方式重复读取同一个元素。

这意味着广播通常能节省内存，但不代表计算一定便宜：

large = np.ones((10_000, 10_000))
bias = np.arange(10_000)

result = large + bias

result 本身仍然需要分配约 800 MB 的 float64 数据空间。广播避免的是 bias 的显式复制，不会消除最终结果的内存需求。

5.5 broadcast_to 返回只读视图

x = np.array([1, 2, 3])
y = np.broadcast_to(x, (2, 3))

print(y)
# [[1 2 3]
#  [1 2 3]]

print(y.flags.writeable)  # False

广播后的多个位置可能对应同一物理元素，因此允许直接写入会产生歧义。若确实需要独立可写数据，应显式复制：

z = np.broadcast_to(x, (2, 3)).copy()
z[0, 0] = 100

5.6 用广播计算两组点之间的距离

设有 m 个二维点和 n 个二维点：

p = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [0.0, 2.0],
])  # (3, 2)

q = np.array([
    [0.0, 1.0],
    [2.0, 2.0],
])  # (2, 2)

希望得到每个 p[i] 与每个 q[j] 的距离，结果形状是 (3, 2)。

先插入轴：

diff = p[:, None, :] - q[None, :, :]
print(diff.shape)  # (3, 2, 2)

形状推导：

p[:, None, :]: (3, 1, 2)
q[None, :, :]: (1, 2, 2)
结果：          (3, 2, 2)

最后一个轴保存二维坐标差：

dist = np.sqrt((diff ** 2).sum(axis=-1))
print(dist.shape)  # (3, 2)

也可以写成：

dist = np.linalg.norm(p[:, None, :] - q[None, :, :], axis=-1)

这种写法避免了显式 Python 双重循环，但中间的 diff 可能很大。若 m 和 n 都很大，需要分块计算，否则广播虽省去了 Python 循环，却可能制造巨大的中间数组。

六、索引：从单个元素到复杂选择

NumPy 索引可以分为：

  1. 基本索引；
  2. 切片；
  3. 布尔索引；
  4. 整数数组索引，也叫高级索引；
  5. 混合索引。

前三类常常产生视图，整数数组索引通常产生副本。不能依赖“看起来像切片”来判断内存共享，应使用 np.shares_memory 验证。

6.1 基本索引和负索引

a = np.array([
    [10, 20, 30],
    [40, 50, 60],
])

print(a[0, 1])  # 20
print(a[-1, -1])  # 60

a[i, j] 表示先沿第 0 轴取第 i 个元素，再沿第 1 轴取第 j 个元素。

print(a[0])      # [10 20 30]
print(a[:, 1])   # [20 50]
print(a[1, :])   # [40 50 60]

: 表示该轴的全部元素。

6.2 切片语法和步长

切片形式为：

start:stop:step

其中 stop 不包含：

x = np.arange(10)

print(x[2:7])    # [2 3 4 5 6]
print(x[::2])    # [0 2 4 6 8]
print(x[::-1])   # [9 8 7 6 5 4 3 2 1 0]

二维切片：

a = np.arange(20).reshape(4, 5)

print(a[1:3, 2:5])
# [[ 7  8  9]
#  [12 13 14]]

切片通常只改变起始地址、形状和步长，不复制元素：

x = np.arange(5)
y = x[1:4]

y[0] = 100
print(x)
# [  0 100   2   3   4]

如果不希望修改原数组：

y = x[1:4].copy()

6.3 视图与副本的根本区别

视图共享底层数据，副本拥有独立数据。

x = np.arange(6)

view = x[::2]
copy = x[::2].copy()

view[0] = 100
copy[1] = 200

print(x)
# [100   1   2   3   4   5]

view[0] 对应原数组的 x[0]，而 copy 与 x 没有共享数据。

可以用以下方式诊断：

print(np.shares_memory(x, view))  # True
print(np.shares_memory(x, copy))  # False

链式索引会增加判断难度：

a[0][1]

通常应优先写成：

a[0, 1]

不仅更清晰，也更容易避免中间对象是视图还是副本的问题。

6.4 布尔索引

布尔索引要求布尔数组能够与被索引区域匹配：

x = np.array([3, 8, 2, 9, 4])
mask = x > 5

print(mask)
# [False  True False  True False]

print(x[mask])
# [8 9]

常见筛选：

selected = x[(x >= 3) & (x <= 8)]

必须使用 &、|、~，并给每个条件加括号：

x[(x > 2) & (x < 9)]

错误写法：

# x > 2 and x < 9

and 和 or 只能处理单个布尔值，不能处理 NumPy 数组。

布尔索引结果通常是新数组：

selected = x[x > 5]
selected[0] = 100

print(x)
# [3 8 2 9 4]

如果需要原地修改满足条件的元素，可以使用赋值：

x[x > 5] = 0
print(x)
# [3 0 2 0 4]

6.5 np.where：条件选择和位置查询

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

result = np.where(x >= 0, x, 0)
print(result)
# [0 0 3 0]

np.where(condition, value_if_true, value_if_false) 会按元素选择结果。

只传一个参数时，它返回满足条件的位置：

indices = np.where(x < 0)
print(indices)
# (array([0, 3]),)

一维数组中常见的等价写法：

np.flatnonzero(x < 0)
# [0 3]

6.6 整数数组索引

整数数组索引按给定位置收集元素：

x = np.array([10, 20, 30, 40, 50])
indices = np.array([4, 0, 2])

print(x[indices])
# [50 10 30]

它不是切片，因此结果通常是副本，并且可以重复索引：

print(x[[1, 1, 3]])
# [20 20 40]

二维数组中：

a = np.array([
    [10, 11, 12],
    [20, 21, 22],
    [30, 31, 32],
])

rows = np.array([0, 2])
cols = np.array([1, 0])

print(a[rows, cols])
# [11 30]

这里是成对选择：

(rows[0], cols[0]) -> a[0, 1] -> 11
(rows[1], cols[1]) -> a[2, 0] -> 30

如果想得到行列组合，需要显式广播：

print(a[rows[:, None], cols])
# [[11 10]
#  [31 30]]

形状为：

rows[:, None]: (2, 1)
cols:          (2)
结果索引形状：  (2, 2)

6.7 take、take_along_axis 和排序索引

对于沿某个轴收集元素，take 更直观：

a = np.arange(12).reshape(3, 4)

print(np.take(a, [0, 2], axis=1))
# [[ 0  2]
#  [ 4  6]
#  [ 8 10]]

当每一行有不同的索引时，使用 take_along_axis：

a = np.array([
    [30, 10, 20],
    [ 3,  1,  2],
])

order = np.argsort(a, axis=1)
sorted_a = np.take_along_axis(a, order, axis=1)

print(sorted_a)
# [[10 20 30]
#  [ 1  2  3]]

argsort 返回排序后的索引，而不是排序后的值。take_along_axis 再按这些索引收集值。

七、向量化：把逐元素逻辑交给数组运算

7.1 什么是向量化

向量化不是简单地“把循环写在一行”，而是把逐元素操作表达为数组运算，由 NumPy 的底层实现执行：

x = np.arange(1_000_000)
y = x * 2 + 1

它与显式 Python 循环的语义等价：

y = np.empty_like(x)

for i in range(x.size):
    y[i] = x[i] * 2 + 1

区别在于，向量化版本把循环从 Python 解释器层移到了底层数组实现。这样通常减少 Python 对象操作、动态类型检查和逐元素函数调用。

但“向量化”不保证所有问题都更快。以下情况可能抵消收益：

  • 创建了多个巨大中间数组；
  • 操作本身很复杂，无法映射为简单 ufunc；
  • 数据规模太小，函数调用开销占主导；
  • 内存带宽而非 CPU 计算成为瓶颈；
  • 使用了对象型数组；
  • 频繁在非连续内存上访问。

7.2 ufunc：通用函数

ufunc 是对数组逐元素执行的函数，例如：

x = np.array([-1.0, 0.0, 1.0, 4.0])

print(np.abs(x))
print(np.sqrt(x + 1))
print(np.exp(x))
print(np.maximum(x, 0))

ufunc 通常支持：

  • 广播；
  • 类型转换；
  • 输出数组复用；
  • 条件执行；
  • 归约。

例如：

x = np.arange(10, dtype=np.float64)
out = np.empty_like(x)

np.multiply(x, 2, out=out)
print(out)

使用 out 可以复用目标缓冲区，减少一次分配：

x = np.arange(5, dtype=np.float64)
np.multiply(x, 2, out=x)
print(x)
# [0. 2. 4. 6. 8.]

原地操作有风险：如果后续仍需要原始数据，就不能覆盖输入。

7.3 条件执行和 where

许多 ufunc 支持 where：

x = np.array([-2.0, 0.0, 4.0])
out = np.full_like(x, np.nan)

np.sqrt(x, out=out, where=x >= 0)
print(out)
# [nan  0.  2.]

这里：

  • out 提供输出缓冲区；
  • where 决定哪些位置执行平方根；
  • 负数位置保留 out 原来的 nan。

如果省略 out，被 where=False 屏蔽的位置可能没有定义好的内容，因此不能把它当作自动填充机制。

7.4 聚合和轴

x = np.array([
    [1, 2, 3],
    [4, 5, 6],
])

print(x.sum())       # 21
print(x.sum(axis=0)) # [5 7 9]
print(x.sum(axis=1)) # [ 6 15]
print(x.mean(axis=0))

常见聚合包括：

x.min()
x.max()
x.mean()
x.std()
x.var()
x.prod()
x.any()
x.all()

空数组和全为 NaN 的数组需要特别处理：

x = np.array([np.nan, np.nan])

print(np.mean(x))     # nan
print(np.nanmean(x))  # 可能产生 RuntimeWarning，并返回 nan

nanmean 会忽略 NaN，但如果没有任何有效值，结果仍然没有实际统计意义。

7.5 einsum：显式表达索引收缩

矩阵乘法：

a = np.array([[1, 2],
              [3, 4]])

b = np.array([[5, 6],
              [7, 8]])

print(np.einsum("ik,kj->ij", a, b))
# [[19 22]
#  [43 50]]

公式是：

Cij=∑kAikBkjC_{ij} = \sum_k A_{ik}B_{kj}

其中：

  • i 对应结果的行；
  • j 对应结果的列；
  • k 是求和轴。

等价写法：

a @ b
np.matmul(a, b)

einsum 的价值在于可以表达更一般的张量收缩、转置和逐元素乘积：

x = np.arange(6).reshape(2, 3)
y = np.ones((2, 3))

row_dot = np.einsum("ij,ij->i", x, y)
print(row_dot)
# [ 3. 12.]

这里：

ri=∑jXijYijr_i = \sum_j X_{ij}Y_{ij}

对于复杂表达式，优先验证形状和索引含义，不要把 einsum 当作缩短代码的技巧。

八、从 Python 循环到数组表达式：一个完整例子

假设有温度数据，形状为：

(samples, features)

每个样本有 3 个特征，需要：

  1. 减去每个特征的均值；
  2. 除以每个特征的标准差；
  3. 将绝对值大于 2 的结果标记为异常。
import numpy as np

data = np.array([
    [10.0, 100.0, 1.0],
    [12.0, 110.0, 2.0],
    [11.0, 105.0, 1.5],
    [50.0, 108.0, 1.2],
])

mean = data.mean(axis=0, keepdims=True)
std = data.std(axis=0, keepdims=True)

normalized = (data - mean) / std
outlier = np.abs(normalized) > 2

print("mean shape:", mean.shape)
print("std shape:", std.shape)
print("normalized:\n", normalized)
print("outlier:\n", outlier)

关键形状：

data:         (4, 3)
mean:         (1, 3)
std:          (1, 3)
normalized:   (4, 3)
outlier:      (4, 3)

因为 mean 和 std 的形状是 (1, 3)，它们会沿样本轴广播到 (4, 3)。如果错误地写成：

mean = data.mean(axis=0)  # shape (3,)

这段计算依然能够工作，因为 (3,) 同样能广播到 (4, 3)。但 keepdims=True 更明确地保留了“每个特征一个统计量”的结构，在推广到三维或更高维数据时也更安全。

若某一列标准差为零：

std[std == 0] = 1

否则会产生除零警告和 NaN 或无穷大。这个修正改变了统计解释：常数列不应被当作普通可标准化特征，工程上也可以选择删除该列或单独处理。

九、内存布局：C-order、Fortran-order 与 strides

9.1 C-order

C-order，也叫行主序，表示最后一个轴变化最快。

a = np.arange(6).reshape(2, 3)

内存逻辑顺序为：

a[0, 0], a[0, 1], a[0, 2],
a[1, 0], a[1, 1], a[1, 2]

对于 int64：

print(a.strides)
# 通常为 (24, 8)

沿最后一个轴移动一格只跳过一个元素，沿第 0 轴移动一格则跳过整行。

9.2 Fortran-order

Fortran-order，也叫列主序，表示第一个轴变化最快：

a = np.asfortranarray(np.arange(6).reshape(2, 3))

print(a.flags.c_contiguous)
print(a.flags.f_contiguous)
print(a.strides)

Fortran-order 的典型步长可能是：

(8, 16)

内存逻辑顺序更接近：

a[0, 0], a[1, 0],
a[0, 1], a[1, 1],
a[0, 2], a[1, 2]

C_CONTIGUOUS 和 F_CONTIGUOUS 表示数组是否符合对应的连续布局：

print(a.flags)

9.3 转置通常改变步长而不是复制数据

x = np.arange(6).reshape(2, 3)
y = x.T

print(x.strides)
print(y.strides)
print(np.shares_memory(x, y))  # True

y 通常只是以相反轴顺序解释同一块数据。因此：

y.flags.c_contiguous

往往为 False。

如果某个底层库要求 C-contiguous 数组，可以显式转换：

y_c = np.ascontiguousarray(y)
print(y_c.flags.c_contiguous)  # True

这可能产生复制，应将它视为有成本的边界操作。

9.4 切片也可能造成非连续访问

x = np.arange(10)
y = x[::2]

print(y)
# [0 2 4 6 8]

print(y.strides)
# 通常为 (16,)；int64 每个元素 8 字节

y 的元素逻辑上连续，但物理上每次跳过一个原数组元素。对这类数组进行大量计算时，内存访问局部性可能比连续数组差。

二维情况：

a = np.arange(16).reshape(4, 4)
b = a[:, ::2]

print(b)
print(b.flags.c_contiguous)

b 只取偶数列，通常不是 C-contiguous。

9.5 内存布局为什么影响性能

现代 CPU 从缓存中按缓存行读取数据。若访问顺序与内存布局一致，连续元素更可能已经被预取到缓存中。

例如 C-order 数组中：

a = np.ones((1000, 1000))

按行访问通常更接近连续内存：

for i in range(a.shape[0]):
    row = a[i, :]
    # 访问 row 中连续元素

按列访问会以较大的步长跳跃：

for j in range(a.shape[1]):
    col = a[:, j]

但这只是常见实现下的性能倾向，不是 NumPy API 对运行时间的规范保证。实际表现还受 BLAS、CPU 缓存、线程、数据类型、数组规模和底层函数实现影响。

十、copy、view 和所有权

10.1 显式创建副本

x = np.arange(5)
y = x.copy()

y[0] = 100
print(x)
# [0 1 2 3 4]

当数据需要跨组件、线程或生命周期边界传递，且生产者可能继续修改原数组时，显式复制可以建立所有权边界。

10.2 view 只改变解释方式

x = np.arange(4, dtype=np.int32)
y = x.view()

y[0] = 100
print(x)
# [100   1   2   3]

view() 不复制底层数据，只创建新的数组对象。

也可以通过不同 dtype 重新解释位模式，但这不是数值转换：

x = np.array([1, 2], dtype=np.int32)
y = x.view(np.uint8)

print(y.shape)
print(y.dtype)

astype 是转换：

z = x.astype(np.float64)

二者语义完全不同：

  • view(dtype)：同一位模式，用另一种方式解释；
  • astype(dtype)：读取原数值并生成目标类型的数据。

10.3 base 只是辅助信息

x = np.arange(5)
y = x[1:]

print(y.base is x)

简单切片中可能如此，但复杂视图不一定直接以原数组作为 base。因此更可靠的判断是：

np.shares_memory(x, y)

如果只需要快速、可能保守地判断，可以使用 np.may_share_memory，但它可能报告潜在共享而非确定共享。

十一、缺失值、掩码和无效数

NumPy 没有统一覆盖所有数据类型的“缺失值”概念。浮点数组通常使用 np.nan：

x = np.array([1.0, np.nan, 3.0])

print(np.isnan(x))
# [False  True False]

print(np.nanmean(x))
# 2.0

整数数组不能直接保存 NaN：

# np.array([1, np.nan], dtype=np.int64)  # 会失败或发生不符合预期的类型处理

可选方案包括：

  • 转成浮点类型；
  • 使用额外布尔掩码；
  • 使用 np.ma.MaskedArray；
  • 在表格数据场景交给 pandas 的缺失值体系。

掩码数组示例：

x = np.ma.array(
    [1.0, 2.0, 999.0],
    mask=[False, False, True],
)

print(x.mean())
# 1.5

999.0 仍然存在于底层数据中，但被掩码标记为无效。掩码不是删除元素，因此需要确认下游函数是否理解掩码语义。

十二、随机数：使用独立生成器

新代码应显式创建随机数生成器：

rng = np.random.default_rng(42)

samples = rng.normal(
    loc=0.0,
    scale=1.0,
    size=(3, 2),
)

print(samples)

这里的 42 是随机种子，用于复现实验。固定种子不代表结果跨所有 NumPy 版本、平台和算法实现都永久一致；它主要用于同一环境下的可重复测试和调试。

常见用法：

rng.integers(0, 10, size=5)
rng.random((2, 3))
rng.choice(["train", "test"], size=10)
rng.permutation(10)

不要在业务函数内部每次调用都重新设置相同种子，否则每次调用可能产生完全相同的随机序列。更好的方式是由调用方创建生成器并传入。

十三、错误处理和诊断

13.1 形状错误

a = np.ones((2, 3))
b = np.ones((2,))

a + b

典型错误是：

ValueError: operands could not be broadcast together

诊断时先打印：

print(a.shape, b.shape)

然后从右向左逐轴比较，而不是盲目插入 reshape。

13.2 索引错误

x = np.arange(3)
x[3]

会产生越界错误：

IndexError

布尔索引的错误通常来自形状不匹配：

x = np.arange(5)
mask = np.array([True, False])

x[mask]

先检查：

print(x.shape)
print(mask.shape)
print(mask.dtype)

13.3 警告不一定是异常

x = np.array([1.0, 0.0])
result = 1 / x

可能产生除零警告，并得到无穷大：

[1. inf]

可以临时控制浮点错误处理：

with np.errstate(divide="raise", invalid="raise"):
    result = 1 / x

此时某些浮点异常会转为 FloatingPointError：

try:
    with np.errstate(divide="raise", invalid="raise"):
        result = 1 / x
except FloatingPointError as exc:
    print("数值计算失败:", exc)

不要在全局范围内永久屏蔽所有警告，否则会隐藏数据质量问题。

13.4 检查有限值

x = np.array([1.0, np.nan, np.inf, -np.inf])

print(np.isfinite(x))
# [ True False False False]

在进入模型、序列化或写入数据库之前，可以显式检查：

if not np.all(np.isfinite(x)):
    raise ValueError("输入包含 NaN 或无穷大")

但是否允许 NaN 要由业务语义决定，不能把所有非有限值都简单视为错误。

十四、性能：真正快在哪里，又会慢在哪里

14.1 Python 循环、数组循环和对象数组

下面的代码在 Python 层逐元素运行：

result = []
for value in values:
    result.append(value * 2 + 1)

数组表达式通常更适合固定类型数值数据：

result = values * 2 + 1

但如果数组是对象型：

objects = np.array([1, 2, 3], dtype=object)

许多操作需要调用 Python 对象的运算协议，无法获得常规数值数组的执行路径。

14.2 避免无意义的中间数组

result = ((x - mean) / std) ** 2

这可能至少产生多个中间结果：

  1. x - mean；
  2. 除以 std；
  3. 平方。

对于中等规模数据，清晰性通常比微优化重要；对于巨大数组，可以考虑复用输出：

result = np.empty_like(x)

np.subtract(x, mean, out=result)
np.divide(result, std, out=result)
np.square(result, out=result)

这种写法减少分配，但增加了原地修改风险，也可能让代码更难审查。只有在内存压力被确认后才值得使用。

14.3 使用基准测试，而不是凭直觉

import timeit

setup = """
import numpy as np
x = np.arange(1_000_000, dtype=np.float64)
"""

vectorized = "y = x * 2 + 1"

print(timeit.timeit(vectorized, setup=setup, number=10))

基准测试应控制：

  • 数据规模；
  • 数据类型；
  • 是否包含数组创建；
  • 是否包含内存分配；
  • 预热和重复次数；
  • 单线程或多线程环境；
  • 是否测量了结果使用。

如果结果数组从未被使用，编译器或运行时环境可能改变实际测量意义；测试应确保计算结果被保留或验证。

十五、与 Python 原生扩展的边界

NumPy 的向量化解决的是“已有操作可以用底层数组循环表达”的问题。遇到无法表达为现有 ufunc 或数组组合的复杂算法时，可以考虑：

  1. 先用纯 Python 验证语义；
  2. 找到数据形状、类型和边界；
  3. 使用 NumPy 表达可向量化部分；
  4. 对剩余热点使用 Cython、PyO3、C API 或其他扩展机制；
  5. 明确 Python 对象与 NumPy 缓冲区之间的所有权和生命周期。

这时最重要的不是“把循环改成 C”，而是定义边界：

Python 对象
    ↓ 参数检查、形状检查、错误转换
NumPy ndarray
    ↓ dtype、shape、strides、flags
底层原生循环
    ↓ 输出缓冲区和异常路径
Python 结果

扩展代码必须明确：

  • 是否要求 C-contiguous；
  • 是否接受只读数组；
  • 是否会保留数组指针；
  • 调用期间 Python 对象是否仍然存活；
  • 是否释放 GIL；
  • 原生异常如何转换为 Python 异常；
  • ABI 是否与目标 Python、NumPy 和编译器兼容。

NumPy 官方用户指南将 C API、F2PY、底层实现和互操作性作为高级主题单独列出。(numpy.org)

一个常见错误是只检查 shape，却忽略 dtype 和 strides：

def accept_matrix(x: np.ndarray) -> None:
    if x.ndim != 2:
        raise ValueError("必须是二维数组")

更完整的边界检查可能是：

def require_float64_c_matrix(x: np.ndarray) -> np.ndarray:
    x = np.asarray(x)

    if x.ndim != 2:
        raise ValueError(f"需要二维数组，实际为 {x.ndim} 维")

    if x.dtype != np.float64:
        x = x.astype(np.float64, copy=False)

    if not x.flags.c_contiguous:
        x = np.ascontiguousarray(x)

    if not x.flags.writeable:
        raise ValueError("底层算法要求可写数组")

    return x

这里 copy=False 只是“不主动要求复制”，并不保证一定不复制；类型不兼容时仍可能生成新数组。ascontiguousarray 在必要时也会复制。

十六、与 pandas 的边界

NumPy 的核心抽象是按位置组织的同质多维数值数组。pandas 的核心抽象是带标签的一维 Series 和二维 DataFrame，强调索引对齐、缺失值、分组、连接和表格语义。pandas 官方文档将其定位为 Python 的数据结构与数据分析工具。(pandas.pydata.org)

例如：

import numpy as np
import pandas as pd

values = np.array([
    [1.0, 10.0],
    [2.0, 20.0],
])

frame = pd.DataFrame(
    values,
    index=["alice", "bob"],
    columns=["age", "score"],
)

NumPy 中：

values[0]

表示按位置取第一行。

pandas 中：

frame.loc["alice"]

表示按标签取名为 "alice" 的行。

如果先把带标签的表格转成 NumPy：

array = frame.to_numpy()

就会丢失行标签和列标签。反向转换时，需要显式传回标签：

pd.DataFrame(
    array,
    index=frame.index,
    columns=frame.columns,
)

因此：

  • 数值内核、矩阵运算、图像张量、科学计算通常适合 NumPy；
  • 需要列名、索引对齐、分组和连接的表格数据通常适合 pandas；
  • 两者转换时必须检查 dtype、缺失值、维度和标签是否仍符合语义。

十七、综合示例：从原始数据到批量计算

下面构造一个小型传感器数据集：

import numpy as np

rng = np.random.default_rng(2026)

# 形状：(样本数, 传感器数)
readings = rng.normal(
    loc=[20.0, 100.0, 0.5],
    scale=[2.0, 10.0, 0.1],
    size=(8, 3),
).astype(np.float64)

print(readings.shape)  # (8, 3)
print(readings.dtype)  # float64

loc 和 scale 的形状都是 (3,)，会广播到 size=(8, 3)，因此每个传感器使用独立的均值和标准差。

计算每个传感器的统计量：

mean = readings.mean(axis=0, keepdims=True)
std = readings.std(axis=0, keepdims=True)

print(mean.shape)  # (1, 3)
print(std.shape)   # (1, 3)

标准化：

safe_std = np.where(std == 0, 1.0, std)
z = (readings - mean) / safe_std

标记异常：

is_outlier = np.abs(z) > 3
outlier_count = is_outlier.sum(axis=0)

print(outlier_count.shape)  # (3,)
print(outlier_count)

计算每个样本的欧氏距离：

distance_from_mean = np.linalg.norm(
    readings - mean,
    axis=1,
)

print(distance_from_mean.shape)  # (8,)

这里：

readings: (8, 3)
mean:      (1, 3)
差值:      (8, 3)
axis=1:    压缩传感器轴
结果:      (8,)

如果要取出异常样本：

sample_is_outlier = is_outlier.any(axis=1)
abnormal_rows = readings[sample_is_outlier]

print(abnormal_rows.shape)

布尔索引会返回副本。如果后续要在原数组上修正异常值，应使用赋值：

readings[is_outlier] = np.nan

但此时数组已经是浮点类型；如果原数组是整数类型，写入 NaN 就不成立。因此缺失值策略必须在创建数组时就考虑。

十八、常见误解与失败模式

误解一：shape 一样就代表数据含义一样

下面两个数组形状相同：

image_hwc.shape == (224, 224, 3)
image_chw.shape == (224, 224, 3)

如果第三个轴在一个数组中表示颜色通道，而在另一个数组中表示宽度，形状虽然一致，语义却完全不同。NumPy 只能检查形状和类型，不能检查业务含义。

误解二：广播失败时随便 reshape

a.shape == (4, 3)
b.shape == (4,)

若意图是给每行加一个值，应写：

a + b[:, None]

这会将 b 变为 (4, 1)。

若意图是给每列加一个值，则 b 必须本来就是 (3,)。强行调整形状可能让程序运行，但结果语义错误。

误解三：切片一定复制，整数索引一定视图

通常情况正好相反：

x[1:4]        # 通常是视图
x[[1, 2, 3]]  # 通常是副本

但复杂组合索引可能改变行为。需要共享内存时检查：

np.shares_memory(original, selected)

误解四：.T 等于“矩阵转置并复制”

.T 通常只是轴顺序变化的视图：

y = x.T

它可能变成非连续数组，也可能导致下游库在边界处复制。若需要独立、连续的数据：

y = x.T.copy(order="C")

误解五：向量化一定比循环快

如果向量化表达式创建了数个巨大中间数组，或者算法本身包含复杂分支，收益可能不明显。对于需要逐元素复杂状态机的逻辑，盲目广播可能造成巨大内存压力。

应先确认：

print(x.shape, x.dtype, x.strides)

再用真实数据进行基准测试，并同时观察内存峰值。

误解六：np.array 总会复制输入

x = np.arange(5)
y = np.array(x)

在某些情况下 NumPy 可以复用已有数组；如果必须保证副本，应写：

y = np.array(x, copy=True)

如果目标是“尽量不复制并转成数组”，常用：

y = np.asarray(x)

但“尽量不复制”不等于“绝不复制”。

十九、建立可靠的 NumPy 心智模型

处理一个数组表达式时，可以按以下顺序分析：

第一步：确认每个对象的语义

print(type(x))
print(x.dtype)
print(x.shape)

不要只看数值内容，要知道每个轴代表什么。

第二步：写出广播前的形状

例如：

data:  (batch, feature)
mean:  (feature,)

再判断它是按特征广播，还是需要：

mean[None, :]

或：

mean[:, None]

第三步：区分轴保留和轴压缩

x.sum(axis=1)
x.sum(axis=1, keepdims=True)

二者数值可能相同，但形状分别是 (n,) 和 (n, 1)，后续广播行为不同。

第四步：确认是否共享内存

view = x[:, ::2]
print(np.shares_memory(x, view))

任何可能修改数据的操作，都应该明确所有权。

第五步：检查连续性和步长

print(x.strides)
print(x.flags.c_contiguous)
print(x.flags.f_contiguous)

当数组要传递给 C、Fortran、GPU 或其他底层库时，dtype、字节序、步长、可写性和连续性都可能成为接口前提。

第六步：评估中间数组

result = f(g(h(x)))

即使最终结果很小，中间结果也可能很大。广播能够避免某些输入复制，但不能自动避免所有计算和输出内存。

二十、结语：从“会调用函数”到“理解数组”

NumPy 的基础不是记住更多函数，而是掌握几个相互连接的机制：

  1. ndarray 由数据、dtype、shape 和 strides 共同描述；
  2. shape 决定轴结构，轴决定聚合、索引和广播方向；
  3. 广播从右向左比较维度，兼容条件是相等或其中一个为 1；
  4. 切片通常创建视图，高级索引通常创建副本；
  5. 向量化把固定类型的逐元素循环交给底层实现，但可能产生中间数组；
  6. C-order、Fortran-order、转置和步长决定数据访问方式；
  7. dtype、可写性和连续性是 NumPy 与原生扩展之间的重要边界；
  8. pandas 在标签和表格语义上提供更高层抽象，转换时不能忽略索引和缺失值语义。

当一个 NumPy 表达式出现错误或性能问题时，最有效的排查入口通常不是继续尝试函数，而是把它还原成：

每个对象的 dtype
每个对象的 shape
每个轴的含义
广播后的目标 shape
是否共享内存
strides 与连续性
中间结果的规模

掌握这些信息后，ndarray、广播、索引、向量化和内存布局不再是互相独立的 API，而会成为同一个数组计算模型的不同表现。

系列导航与关联阅读

官方资料

本文依据 Python 官方文档、相关 PEP 与生态项目官方文档重新梳理；正文、示例与工程清单由 WR BLOG 编写。