NumPy这个库上篇我们聊完了它的热身动作——怎么创建ndarray、切片、reshape、拼接拆分这些都是“用手搬砖”的阶段。但真正让我从“会用NumPy”变成“觉得NumPy真香”的是下篇这些内容广播机制、高级索引、通用函数、线性代数、随机数、结构化数组以及视图和副本背后那些隐藏的内存逻辑。这一篇笔记就是这些内容的完整扩展。先说我为什么写这篇。最近在做一个多站点气象数据的预处理数据形状大概是(时间点, 站点数, 高度层)。需求很朴素每个站点、每个高度层分别做z-score标准化也就是 (x - mean) / std。第一版我老老实实写三层for循环跑了十几秒不说代码嵌套得自己都不想看。后来换成广播ufunc的写法三行搞定速度直接快了两个数量级。这种体验才是NumPy真正值钱的地方。这篇笔记定位是NumPy的“进阶实战”适合已经熟悉ndarray基本用法、但还没系统接触过广播、ufunc、性能调优的Python学习者。我会把每个主题都落到一个可以直接跑的小例子上还会穿插一些我踩过的坑和简单的性能对比。内容偏底层但绝不枯燥因为每个知识点都有对应的实战场景。1. 广播机制让不同shape的数组直接做运算的规则1.1 从一次shape报错说起广播的底层规则写代码时最常看见的报错之一就是“operands could not be broadcast together with shapes ...”。两个数组相加shape对不上NumPy就报这个错。理解广播就是理解NumPy什么时候允许“对不上”以及它怎么处理这种“对不上”。广播的规则其实只有一句话从最右边的维度开始往左对齐依次比较每个维度的大小满足两条之一就可以运算两个维度相等其中一个维度是1其中一个数组在这个维度上根本没有也就是维度缺失相当于隐式补1。比较的维度只要出现一次“不满足上面任何一个条件”整个运算就报错。提示这里有个容易混淆的地方维度缺失是在“最左边”补1所以(4,)会自动变成(1, 4)但不会变成(4, 1)。想要(4, 1)必须手动加一个长度为1的轴。这个规则很多人背得下来但用的时候还是会懵。我自己的记忆方式是把它类比成一个“自动补小纸片”的过程小的数组会在缺少维度的方向上自动扩展扩展的那一格不是真复制数据而是重复用同一个值。比如一个形状为(3, 1)的数组和一个形状为(1, 4)的数组相加NumPy会把前者的“1”那一维扩展成4把后者的“1”那一维扩展成3最终结果形状是(3, 4)。import numpy as np a np.array([[1], [2], [3]]) # 形状 (3, 1) b np.array([[10, 20, 30, 40]]) # 形状 (1, 4) c a b print(c.shape) # (3, 4) print(c)输出(3, 4) [[11 21 31 41] [12 22 32 42] [13 23 33 43]]可以看到a的每一行被横向复制了4次b的每一行被纵向复制了3次然后逐元素相加。关键是这个“复制”是逻辑上的内存里并没有真的复制出12个元素所以广播在空间上非常省。1.2 三种高频广播场景和对应的内存行为把广播用熟之后你会发现日常80%的数组运算都可以归到三类场景里数组与标量运算data * 2、data - mean、data / std。标量被认为是形状为()的数组广播到和data完全相同的形状。沿某个轴减去一个统计量最典型的例子是行方向的均值中心化。假设X的形状是(样本数, 特征数)你想对每个特征减去该特征的均值代码是X - X.mean(axis0)。这里的X.mean(axis0)形状是(特征数,)也就是(4,)NumPy会自动把它当作(1, 4)然后逐行广播到(样本数, 特征数)。网格坐标生成np.arange(5)[:, None] * np.arange(4)[None, :]通过引入长度为1的维度生成一个5x4的网格结果。这种写法在画等高线、计算距离矩阵时特别常见。这里有一个容易忽略的点X - X.mean(axis0)这种写法之所以能工作是因为一维数组(4,)在广播时默认被放在“最右边”即当作行向量。如果你需要的是把一维数组(4,)当列向量用就必须显式改成X - X.mean(axis1)[:, None]否则NumPy会在维度上对不齐结果完全不是你想要的那个。1.3 广播实战案例行归一化与距离矩阵的高效写法举一个我在项目中真正用到的例子。一批样本数据X形状(1000, 50)要做标准化每个特征减均值除标准差。用广播可以写成mean X.mean(axis0) # (50,) std X.std(axis0) # (50,) X_norm (X - mean) / std同样的逻辑如果写成for循环需要两重循环代码长且容易出错。广播版本一次到位而且因为底层是C循环速度上差距不是一点半点。另一个例子是计算样本间的欧氏距离矩阵。假设有两个点集A和BA形状(m, d)B形状(n, d)要得到距离矩阵D形状(m, n)其中D[i, j] ||A[i] - B[j]||。用广播加np.sum可以这样写D np.sqrt(((A[:, None, :] - B[None, :, :]) ** 2).sum(axis-1))A[:, None, :]把A扩展成(m, 1, d)B[None, :, :]把B扩展成(1, n, d)两者相减时就会得到一个(m, n, d)的差张量每个位置对应一对点之间的坐标差向量。再对最后一维求平方和开根号距离矩阵就出来了。这段代码第一次看会觉得“玄乎”但只要掌握了“在中间插一个长度为1的维度让广播自动对齐”这个思路就很容易举一反三。后面讲性能优化时我还会再用这个例子说明向量化和Python循环的速度差异。2. 高级索引与掩码筛选把数据分析的活干利索2.1 布尔掩码索引条件筛选的正解NumPy里最常用、也最直接的高级索引是布尔掩码。arr[arr 5]这一句会返回一个一维数组里面是所有大于5的元素。这里的arr 5本身就是一个布尔数组把它当作索引传给原数组NumPy就会挑出所有对应位置为True的元素。arr np.array([12, 3, 8, 15, 6, 9]) mask (arr 5) (arr 12) print(arr[mask]) # [ 8 6 9]这里有个细节值得特别提醒布尔掩码可以和逻辑运算符自由组合但必须用与、|或、~非而不是Python的and/or/not。原因很朴素and/or要求操作数是单个布尔值而arr 5是一个数组Python的解释器会尝试调用它的__bool__方法数组的__bool__不明确于是直接报“truth value of an array is ambiguous”。另外的优先级比低所以保险起见每个条件外面都要加括号。我见过不少新手因为漏了括号得到的mask永远是一个由逐位运算得到的“意外结果”排查半天才发现是优先级问题。布尔掩码的一个漂亮特性是它天然适合用来做数据清洗。比如有一个数组里混入了无效值-9999想全部去掉只需valid data[data ! -9999]一行搞定。2.2 花式索引的取值逻辑与副本陷阱花式索引指的是把一个整数数组作为索引arr[[0, 2, 2, -1]]。它和普通切片最大的区别是它会返回一个新数组而不是原数组的视图。而且索引可以重复可以乱序结果数组的形状由索引数组的形状决定而不是原数组的形状。arr np.arange(10) print(arr[[0, 2, 2, -1]]) # [0 2 2 9]对于二维数组花式索引可以同时指定行索引和列索引arr[[0, 1, 3], [2, 2, 0]]会返回三个元素(0,2)、(1,2)、(3,0)并不是返回一个子矩阵。如果你想要子矩阵需要写成arr[np.ix_([0, 1, 3], [2, 2, 0])]或者用两次花式索引。np.ix_这个工具很少被提到但在需要取“某些行和某些列的交叉块”时非常实用。需要注意一个坑花式索引返回的是副本。也就是说arr[[0, 1, 2]] something这样写是能改原数组的因为赋值操作是通过__setitem__实现的但如果先写sub arr[[0, 1, 2]]然后去改sub原数组不会有任何变化。在debug半天怀疑“我明明改了为什么数组没变”的时候多半就是栽在这个副本与视图的区别上。2.3 np.where与np.select按条件重编码的高效实现np.where(condition, x, y)是一个向量化的三目运算符condition为True的位置取xFalse的位置取y。三个参数都不传x和y时np.where(condition)等价于np.nonzero(condition)返回满足条件的下标元组这个功能也常用来定位异常值。scores np.array([45, 78, 92, 58, 66]) label np.where(scores 60, pass, fail) print(label) # [fail pass pass fail pass]当分支超过两个时np.where可以嵌套但嵌套可读性很差。更好的方案是np.select(condlist, choicelist)两个参数都是列表按顺序匹配第一个满足的条件生效都不满足时用默认值condlist [scores 60, scores 85] choicelist [不及格, 良好] label np.select(condlist, choicelist, default优秀) print(label) # [不及格 良好 优秀 不及格 良好]这段代码比嵌套的np.where清晰得多也更容易维护。实际业务里常见的“把连续变量分箱成类别变量”这类需求我基本都用np.select来做。关于性能如果条件判断的次数不多、数据量也不大直接写Python的列表推导也没问题。但数据量一旦到百万级别以上np.where和np.select的优势就很明显因为它们把循环下沉到了C层。3. 通用函数ufunc把循环写进C层3.1 ufunc的本质为什么比Python循环快那么多ufunc的全称是universal function翻译成“通用函数”但它更贴近的理解是“对ndarray逐元素运算的函数”。加、减、乘、除、幂、指数、对数、三角函数、比较运算这些都是ufunc。np.add对应np.multiply对应*np.greater对应。平时写arr * 2其实就是在调用np.multiply(arr, 2)。ufunc为什么快因为它避免了Python解释器逐元素循环的开销。Python的for循环每一步都要做类型检查、解释执行而ufunc直接在C层面遍历内存中连续排列的数据一次循环完成整个数组的运算。这也是为什么“能用向量化就不要用for循环”成为NumPy社区的铁律。a np.arange(1_000_000) # 推荐 b np.sin(a) np.cos(a) # 不推荐 b np.array([np.sin(x) np.cos(x) for x in a])后者的耗时通常是前者的几十到上百倍。我不是说Python循环一无是处——在数据量小、逻辑复杂、不可向量化的时候循环反而是清晰的。但凡是能用ufunc表达的逻辑优先用ufunc这是写NumPy代码的基本素养。3.2 reduce、accumulate、outerufunc的隐藏接口ufunc除了直接逐元素运算还带了一些聚合方法这些方法平时被隐藏得很深但用起来非常惊艳。ufunc.reduce把整个数组按某个轴聚合成一个值。np.add.reduce(arr)就是arr.sum()np.maximum.reduce(arr)就是arr.max(axis0)。在逻辑运算时它表现得非常自然np.logical_and.reduce(bool_array)可以一次性判断“所有位置是否都为True”。ufunc.accumulate累积运算。np.add.accumulate(arr)返回前缀和每个位置的值是原数组从开始到当前位置的和。这个操作在计算累计收益、滑动统计量时很方便。ufunc.outer外积。np.multiply.outer(a, b)生成a[i] * b[j]的矩阵。np.add.outer(a, b)则是一个“两两相加”的矩阵在动态规划、网格搜索里很常见。returns np.array([0.01, -0.02, 0.03, 0.005, -0.01]) nav np.add.accumulate(np.log1p(returns)) # 对数累计 nav np.exp(nav) print(nav) # [1.01005 0.99005 1.02007 1.02517 1.01493]一行代码画出从1开始的累计净值曲线不用写for循环去逐个累乘。类似这种“序列型”操作accumulate几乎总是最优解。3.3 自定义ufuncnp.frompyfunc与np.vectorize的边界有时候标准ufunc覆盖不了我们的业务逻辑比如要处理一个带复杂规则的函数。这时可以用np.frompyfunc把一个Python函数包装成ufunc风格的函数也可以用它来给函数指定输入输出参数个数。def classify(x): if x 0: return pos elif x 0: return neg return zero ufunc_classify np.frompyfunc(classify, 1, 1) arr np.array([-1, 0, 2, 3, -5]) print(ufunc_classify(arr)) # [neg zero pos pos neg]np.vectorize也提供类似功能但它本质上是包装了一层循环并不是真正的C加速。所以我的建议是如果你确实需要自定义逐元素函数先检查能不能用np.where、np.select或者标准ufunc的组合来表达实在不行再用frompyfunc或vectorize但要有心理准备它的速度不会比Python循环快太多只是代码简洁一些。顺带一提自定义ufunc返回的对象类型通常是object需要数值运算时记得转一下astype。4. 线性代数与随机数科学计算的左右手4.1 np.linalg模块速览先记住这几个函数NumPy的np.linalg模块提供了常见线性代数运算。先列个速查表函数作用典型场景np.linalg.inv矩阵求逆解线性系统、统计计算np.linalg.det行列式判断矩阵是否奇异、体积计算np.linalg.solve解线性方程组比先求逆再乘更稳定且更快np.linalg.eig特征值和特征向量PCA主成分分析np.linalg.svd奇异值分解降维、推荐系统、图像压缩np.linalg.norm向量/矩阵范数距离计算、正则化np.linalg.qrQR分解最小二乘、数值稳定的求特征值这里最值得强调的一点是解线性方程组Ax b时优先用np.linalg.solve(A, b)而不是np.linalg.inv(A) b。原因倒不是慢而是inv先显式计算了逆矩阵再乘b中间会引入额外的浮点误差而且计算量更大。solve内部用的是LU分解数值稳定性更好速度也更快。4.2 行列式与矩阵求逆的数值稳定性一个容易被忽略的坑np.linalg.det用来求行列式。行列式接近0意味着矩阵接近奇异也就是矩阵的某些行或列近乎线性相关这种情况下求逆的结果会在数值上非常不稳定。实际业务中行列式很难严格等于0往往是一个很小的值这时候inv的结果会出现巨大的数值波动。我之前处理过一个最小二乘问题设计矩阵的条件数特别大直接求逆得到的结果和用np.linalg.lstsq得到的结果差了三个数量级。查资料才发现线性最小二乘的“教科书公式”X^T X的逆再乘X^T在数值上是不稳定的。正确的做法是优先用np.linalg.lstsq或者np.linalg.pinv伪逆它们在内部会做SVD分解对接近奇异的情况处理得更好。所以给大家一个实操建议凡是涉及“求逆”的场景先想一想能不能换成solve、lstsq、pinv中的某一个。大多数情况下求逆只是数学推导里的中间步骤落到代码上并不需要真的算逆矩阵。4.3 np.random的进化从RandomState到Generator随机数生成也是NumPy中日常频率很高的部分。老版本的写法是np.random.seed(42)加np.random.randn(...)、np.random.randint(...)这套基于RandomState的接口虽然短但存在一些问题全局种子容易干扰其他模块的随机性而且随机数算法的选择不够灵活。NumPy 1.17开始推荐的新接口是Generator用法如下rng np.random.default_rng(42) # 自己的随机数发生器 x rng.normal(0, 1, size(3, 3)) # 标准正态分布 y rng.integers(0, 100, size10) # 整数随机数 idx rng.choice(100, size5, replaceFalse) # 不放回抽样Generator相比RandomState的几个好处随机数生成速度更快支持PCG64等更新更强的算法用起来更符合“每个过程一个随机源”的工程实践避免全局状态污染。如果你管理多个数据管线建议给每个阶段都初始化一个独立的default_rng而不是共用一个种子。随机数在实际中还能做很多有趣的事。比如用蒙特卡洛方法估算圆周率在一个边长为2的正方形内随机撒点统计落在内切圆内的比例乘以4就是圆周率的近似值rng np.random.default_rng(42) n 100_000 x rng.uniform(-1, 1, n) y rng.uniform(-1, 1, n) pi_approx 4 * ((x**2 y**2) 1).mean() print(pi_approx) # 约 3.14这段代码充分体现了“向量化随机数”组合的威力一个for循环都没有。5. 结构化数组与文件IO让ndarray带上语义5.1 结构化数组给每列起名字的ndarrayNumPy的ndarray默认是同构的也就是每个元素类型一致。但有一种情况例外结构化数组。它允许一个数组里包含多个字段每个字段有自己的名字和类型。可以把它理解成“更简洁的DataFrame雏形”。定义结构化数组非常简单只要在dtype里指定字段名和类型data np.array( [(1, 2.5, bapple), (2, 3.6, bbanana)], dtype[(id, i4), (score, f8), (name, S10)] ) print(data[score]) # [2.5 3.6] print(data[0][name]) # bapple有了字段名之后取整列数据就变成了data[score]不再需要记住“第2列是score”。这种代码在可读性上的提升非常大尤其是当列数多、字段语义重要的时候。结构化数组也可以和普通ndarray互相转换。如果你有一堆csv数据每列的语义各不相同用结构化数组存储会比用一个大的object数组清晰得多。不过要注意结构化数组的算力支持不如普通ndarray全面——很多ufunc是逐元素作用于整个结构体而不是按字段分别作用所以遇到复杂运算时通常还是建议把结构化数组拆成多个普通数组算完再合回结构化数组。5.2 np.loadtxt与np.genfromtxt文本数据读写的细节参数读写文本文件最常见的是np.loadtxt和np.savetxt。np.savetxt保存二维数组np.loadtxt读回来中间用fmt控制格式、delimiter控制分隔符arr np.array([[1.5, 2.5], [3.5, 4.5]]) np.savetxt(data.csv, arr, delimiter,, fmt%.6f) arr2 np.loadtxt(data.csv, delimiter,)loadtxt还有几个参数容易被忽略但很实用usecols只读某些列skiprows跳过文件头unpackTrue把结果按列拆开相当于转置输出。如果数据文件带有列名表头想直接按名字取列更合适的工具是np.genfromtxtarr np.genfromtxt(data_with_header.csv, delimiter,, namesTrue) print(arr[temperature])namesTrue表示第一行是列名读出来的结构化数组可以直接用字段名访问。genfromtxt还有一个好处支持处理缺失值missing_values和filling_values参数可以在读取阶段就把无效数据替换成指定值。5.3 np.savez与np.memmap处理超大数据集的两把钥匙当中间结果太多想一次性存下多个数组用np.savez或np.savez_compressed比逐个np.save方便得多np.savez_compressed(result.npz, XX, yy, labelslabels) loaded np.load(result.npz) X2 loaded[X].npz文件本质是一个zip包savez_compressed会做压缩对重复值多的数据效果尤其明显。再往外走一步如果数据大到根本装不进内存np.memmap就派上用场了。memmap把磁盘上的二进制文件映射到内存地址空间读写时操作系统按页加载不需要一次性读入整个数组。这在处理超大矩阵时非常关键# 创建一个 10000x10000 的 float64 文件约 800MB mm np.memmap(large.dat, dtypefloat64, modew, shape(10000, 10000)) # 可以像普通数组一样按块操作 mm[:, 0] 1 mm.flush()memmap的坑在于它默认返回的是一个memory-mapped数组很多操作会触发“物化”把整块数据读入内存比如np.array(mm)就会真正分配内存。所以用它时要有意识地做分块计算避免无意中把整个数组load进RAM。6. 视图vs副本与性能优化NumPy代码提速的几个实操点6.1 视图还是副本一个切片引发的“血案”对ndarray做基本切片arr[1:5]时返回的是原数组的视图。视图和数据是共享内存的所以修改视图会同步改动原数组。这对性能是好事避免了不必要的复制但对不熟悉的人来说这是一个经典的“血案”来源。a np.arange(10) b a[2:8] b[0] 999 print(a) # [ 0 1 999 3 4 5 6 7 8 9]上面这段代码里我只改了b结果a也跟着变了。很多新手写代码时以为切片是“复制了一份”结果无意间改了原数据排查半天。那么什么操作返回视图什么操作返回副本呢操作返回类型说明基础切片arr[1:5]视图与原数组共享内存整数/布尔索引arr[[0,2]]、arr[mask]副本新分配内存arr.reshape()视图可能在原数据连续时可共享内存arr.ravel()视图优先尽量返回视图arr.flatten()副本总是复制arr.copy()副本显式复制判断两个数组是否共享内存可以用np.shares_memory(a, b)。如果你确实需要一份独立的副本请显式调用.copy()这是最稳妥的做法。6.2 C-order与F-order为什么遍历顺序会影响性能NumPy数组在内存中的排列顺序有行优先C-order和列优先F-order之分。默认是C-order即先沿着行方向填满再换下一行F-order则先填满列方向。这个顺序对性能的影响在于CPU缓存和分页的局部性。如果你的算法是逐行遍历数组C-order会非常顺畅因为相邻元素在内存中也是相邻的如果算法是逐列遍历F-order会更快。用np.ascontiguousarray(arr)可以强制转换成C-ordernp.asfortranarray(arr)转成F-order。举一个典型的例子在按列累加时C-order的数组访问跨越大步长缓存命中率下降性能可能差好几倍。之前我做矩阵乘法优化时把其中一个矩阵转成F-order后耗时下降了约30%。所以如果你的算法有明确的遍历方向检查一下数组的arr.flags[C_CONTIGUOUS]或arr.flags[F_CONTIGUOUS]是值得养成习惯的一件事。6.3 实测提速案例用广播重构三层for循环最后用一个真实案例来收尾。假设我要计算两个点集A和B之间的平方距离矩阵A形状(500, 3)B形状(800, 3)。第一版用Python三层循环写严格说是两层循环因为每对点内部还要遍历3个坐标D np.empty((500, 800)) for i in range(500): for j in range(800): s 0 for k in range(3): s (A[i, k] - B[j, k]) ** 2 D[i, j] s第二版用广播加ufuncdiff A[:, None, :] - B[None, :, :] # (500, 800, 3) D (diff ** 2).sum(axis-1)实测下来前者在我的机器上大约耗时0.8秒左右后者只有约3毫秒速度提升了两百多倍。这个差距不是“稍微优化”的程度而是“明明可以秒开却要卡很久”的差别。还有一个中间的写法用scipy.spatial.distance.cdist它在底层调了向量化的C实现速度也很快代码更简洁。如果项目里已经装了scipy优先用专业工具也是完全合理的选择。不过NumPy版本的优点是不依赖额外库逻辑也足够直观适合作为理解距离计算的模板。最后再说一个我自己的习惯。写NumPy代码时我会先问自己三个问题这个操作能不能用广播写能不能用ufunc写是不是在无意中创建了副本带着这三个问题去看代码很多“慢”都能在没开始跑之前就规避掉。另外如果遇到UserWarning: failed to initialize numpy这类报错多半是环境里存在多个NumPy版本冲突或者某个第三方库引用了不同路径的NumPy。处理方式也很简单在干净环境里重装一遍NumPy或者用pip install --upgrade numpy统一版本。别慌这类问题大多不是代码逻辑的问题而是环境管理的问题。希望这篇笔记能帮你在NumPy上少走点弯路把省下来的时间留给真正需要思考的地方。 
长河君 · 2026/10/3 1:09:22 · 阅读 1.2k NumPy进阶实战：广播机制、向量化与性能优化全解析
NumPy这个库上篇我们聊完了它的热身动作——怎么创建ndarray、切片、reshape、拼接拆分这些都是“用手搬砖”的阶段。但真正让我从“会用NumPy”变成“觉得NumPy真香”的是下篇这些内容广播机制、高级索引、通用函数、线性代数、随机数、结构化数组以及视图和副本背后那些隐藏的内存逻辑。这一篇笔记就是这些内容的完整扩展。先说我为什么写这篇。最近在做一个多站点气象数据的预处理数据形状大概是(时间点, 站点数, 高度层)。需求很朴素每个站点、每个高度层分别做z-score标准化也就是 (x - mean) / std。第一版我老老实实写三层for循环跑了十几秒不说代码嵌套得自己都不想看。后来换成广播ufunc的写法三行搞定速度直接快了两个数量级。这种体验才是NumPy真正值钱的地方。这篇笔记定位是NumPy的“进阶实战”适合已经熟悉ndarray基本用法、但还没系统接触过广播、ufunc、性能调优的Python学习者。我会把每个主题都落到一个可以直接跑的小例子上还会穿插一些我踩过的坑和简单的性能对比。内容偏底层但绝不枯燥因为每个知识点都有对应的实战场景。1. 广播机制让不同shape的数组直接做运算的规则1.1 从一次shape报错说起广播的底层规则写代码时最常看见的报错之一就是“operands could not be broadcast together with shapes ...”。两个数组相加shape对不上NumPy就报这个错。理解广播就是理解NumPy什么时候允许“对不上”以及它怎么处理这种“对不上”。广播的规则其实只有一句话从最右边的维度开始往左对齐依次比较每个维度的大小满足两条之一就可以运算两个维度相等其中一个维度是1其中一个数组在这个维度上根本没有也就是维度缺失相当于隐式补1。比较的维度只要出现一次“不满足上面任何一个条件”整个运算就报错。提示这里有个容易混淆的地方维度缺失是在“最左边”补1所以(4,)会自动变成(1, 4)但不会变成(4, 1)。想要(4, 1)必须手动加一个长度为1的轴。这个规则很多人背得下来但用的时候还是会懵。我自己的记忆方式是把它类比成一个“自动补小纸片”的过程小的数组会在缺少维度的方向上自动扩展扩展的那一格不是真复制数据而是重复用同一个值。比如一个形状为(3, 1)的数组和一个形状为(1, 4)的数组相加NumPy会把前者的“1”那一维扩展成4把后者的“1”那一维扩展成3最终结果形状是(3, 4)。import numpy as np a np.array([[1], [2], [3]]) # 形状 (3, 1) b np.array([[10, 20, 30, 40]]) # 形状 (1, 4) c a b print(c.shape) # (3, 4) print(c)输出(3, 4) [[11 21 31 41] [12 22 32 42] [13 23 33 43]]可以看到a的每一行被横向复制了4次b的每一行被纵向复制了3次然后逐元素相加。关键是这个“复制”是逻辑上的内存里并没有真的复制出12个元素所以广播在空间上非常省。1.2 三种高频广播场景和对应的内存行为把广播用熟之后你会发现日常80%的数组运算都可以归到三类场景里数组与标量运算data * 2、data - mean、data / std。标量被认为是形状为()的数组广播到和data完全相同的形状。沿某个轴减去一个统计量最典型的例子是行方向的均值中心化。假设X的形状是(样本数, 特征数)你想对每个特征减去该特征的均值代码是X - X.mean(axis0)。这里的X.mean(axis0)形状是(特征数,)也就是(4,)NumPy会自动把它当作(1, 4)然后逐行广播到(样本数, 特征数)。网格坐标生成np.arange(5)[:, None] * np.arange(4)[None, :]通过引入长度为1的维度生成一个5x4的网格结果。这种写法在画等高线、计算距离矩阵时特别常见。这里有一个容易忽略的点X - X.mean(axis0)这种写法之所以能工作是因为一维数组(4,)在广播时默认被放在“最右边”即当作行向量。如果你需要的是把一维数组(4,)当列向量用就必须显式改成X - X.mean(axis1)[:, None]否则NumPy会在维度上对不齐结果完全不是你想要的那个。1.3 广播实战案例行归一化与距离矩阵的高效写法举一个我在项目中真正用到的例子。一批样本数据X形状(1000, 50)要做标准化每个特征减均值除标准差。用广播可以写成mean X.mean(axis0) # (50,) std X.std(axis0) # (50,) X_norm (X - mean) / std同样的逻辑如果写成for循环需要两重循环代码长且容易出错。广播版本一次到位而且因为底层是C循环速度上差距不是一点半点。另一个例子是计算样本间的欧氏距离矩阵。假设有两个点集A和BA形状(m, d)B形状(n, d)要得到距离矩阵D形状(m, n)其中D[i, j] ||A[i] - B[j]||。用广播加np.sum可以这样写D np.sqrt(((A[:, None, :] - B[None, :, :]) ** 2).sum(axis-1))A[:, None, :]把A扩展成(m, 1, d)B[None, :, :]把B扩展成(1, n, d)两者相减时就会得到一个(m, n, d)的差张量每个位置对应一对点之间的坐标差向量。再对最后一维求平方和开根号距离矩阵就出来了。这段代码第一次看会觉得“玄乎”但只要掌握了“在中间插一个长度为1的维度让广播自动对齐”这个思路就很容易举一反三。后面讲性能优化时我还会再用这个例子说明向量化和Python循环的速度差异。2. 高级索引与掩码筛选把数据分析的活干利索2.1 布尔掩码索引条件筛选的正解NumPy里最常用、也最直接的高级索引是布尔掩码。arr[arr 5]这一句会返回一个一维数组里面是所有大于5的元素。这里的arr 5本身就是一个布尔数组把它当作索引传给原数组NumPy就会挑出所有对应位置为True的元素。arr np.array([12, 3, 8, 15, 6, 9]) mask (arr 5) (arr 12) print(arr[mask]) # [ 8 6 9]这里有个细节值得特别提醒布尔掩码可以和逻辑运算符自由组合但必须用与、|或、~非而不是Python的and/or/not。原因很朴素and/or要求操作数是单个布尔值而arr 5是一个数组Python的解释器会尝试调用它的__bool__方法数组的__bool__不明确于是直接报“truth value of an array is ambiguous”。另外的优先级比低所以保险起见每个条件外面都要加括号。我见过不少新手因为漏了括号得到的mask永远是一个由逐位运算得到的“意外结果”排查半天才发现是优先级问题。布尔掩码的一个漂亮特性是它天然适合用来做数据清洗。比如有一个数组里混入了无效值-9999想全部去掉只需valid data[data ! -9999]一行搞定。2.2 花式索引的取值逻辑与副本陷阱花式索引指的是把一个整数数组作为索引arr[[0, 2, 2, -1]]。它和普通切片最大的区别是它会返回一个新数组而不是原数组的视图。而且索引可以重复可以乱序结果数组的形状由索引数组的形状决定而不是原数组的形状。arr np.arange(10) print(arr[[0, 2, 2, -1]]) # [0 2 2 9]对于二维数组花式索引可以同时指定行索引和列索引arr[[0, 1, 3], [2, 2, 0]]会返回三个元素(0,2)、(1,2)、(3,0)并不是返回一个子矩阵。如果你想要子矩阵需要写成arr[np.ix_([0, 1, 3], [2, 2, 0])]或者用两次花式索引。np.ix_这个工具很少被提到但在需要取“某些行和某些列的交叉块”时非常实用。需要注意一个坑花式索引返回的是副本。也就是说arr[[0, 1, 2]] something这样写是能改原数组的因为赋值操作是通过__setitem__实现的但如果先写sub arr[[0, 1, 2]]然后去改sub原数组不会有任何变化。在debug半天怀疑“我明明改了为什么数组没变”的时候多半就是栽在这个副本与视图的区别上。2.3 np.where与np.select按条件重编码的高效实现np.where(condition, x, y)是一个向量化的三目运算符condition为True的位置取xFalse的位置取y。三个参数都不传x和y时np.where(condition)等价于np.nonzero(condition)返回满足条件的下标元组这个功能也常用来定位异常值。scores np.array([45, 78, 92, 58, 66]) label np.where(scores 60, pass, fail) print(label) # [fail pass pass fail pass]当分支超过两个时np.where可以嵌套但嵌套可读性很差。更好的方案是np.select(condlist, choicelist)两个参数都是列表按顺序匹配第一个满足的条件生效都不满足时用默认值condlist [scores 60, scores 85] choicelist [不及格, 良好] label np.select(condlist, choicelist, default优秀) print(label) # [不及格 良好 优秀 不及格 良好]这段代码比嵌套的np.where清晰得多也更容易维护。实际业务里常见的“把连续变量分箱成类别变量”这类需求我基本都用np.select来做。关于性能如果条件判断的次数不多、数据量也不大直接写Python的列表推导也没问题。但数据量一旦到百万级别以上np.where和np.select的优势就很明显因为它们把循环下沉到了C层。3. 通用函数ufunc把循环写进C层3.1 ufunc的本质为什么比Python循环快那么多ufunc的全称是universal function翻译成“通用函数”但它更贴近的理解是“对ndarray逐元素运算的函数”。加、减、乘、除、幂、指数、对数、三角函数、比较运算这些都是ufunc。np.add对应np.multiply对应*np.greater对应。平时写arr * 2其实就是在调用np.multiply(arr, 2)。ufunc为什么快因为它避免了Python解释器逐元素循环的开销。Python的for循环每一步都要做类型检查、解释执行而ufunc直接在C层面遍历内存中连续排列的数据一次循环完成整个数组的运算。这也是为什么“能用向量化就不要用for循环”成为NumPy社区的铁律。a np.arange(1_000_000) # 推荐 b np.sin(a) np.cos(a) # 不推荐 b np.array([np.sin(x) np.cos(x) for x in a])后者的耗时通常是前者的几十到上百倍。我不是说Python循环一无是处——在数据量小、逻辑复杂、不可向量化的时候循环反而是清晰的。但凡是能用ufunc表达的逻辑优先用ufunc这是写NumPy代码的基本素养。3.2 reduce、accumulate、outerufunc的隐藏接口ufunc除了直接逐元素运算还带了一些聚合方法这些方法平时被隐藏得很深但用起来非常惊艳。ufunc.reduce把整个数组按某个轴聚合成一个值。np.add.reduce(arr)就是arr.sum()np.maximum.reduce(arr)就是arr.max(axis0)。在逻辑运算时它表现得非常自然np.logical_and.reduce(bool_array)可以一次性判断“所有位置是否都为True”。ufunc.accumulate累积运算。np.add.accumulate(arr)返回前缀和每个位置的值是原数组从开始到当前位置的和。这个操作在计算累计收益、滑动统计量时很方便。ufunc.outer外积。np.multiply.outer(a, b)生成a[i] * b[j]的矩阵。np.add.outer(a, b)则是一个“两两相加”的矩阵在动态规划、网格搜索里很常见。returns np.array([0.01, -0.02, 0.03, 0.005, -0.01]) nav np.add.accumulate(np.log1p(returns)) # 对数累计 nav np.exp(nav) print(nav) # [1.01005 0.99005 1.02007 1.02517 1.01493]一行代码画出从1开始的累计净值曲线不用写for循环去逐个累乘。类似这种“序列型”操作accumulate几乎总是最优解。3.3 自定义ufuncnp.frompyfunc与np.vectorize的边界有时候标准ufunc覆盖不了我们的业务逻辑比如要处理一个带复杂规则的函数。这时可以用np.frompyfunc把一个Python函数包装成ufunc风格的函数也可以用它来给函数指定输入输出参数个数。def classify(x): if x 0: return pos elif x 0: return neg return zero ufunc_classify np.frompyfunc(classify, 1, 1) arr np.array([-1, 0, 2, 3, -5]) print(ufunc_classify(arr)) # [neg zero pos pos neg]np.vectorize也提供类似功能但它本质上是包装了一层循环并不是真正的C加速。所以我的建议是如果你确实需要自定义逐元素函数先检查能不能用np.where、np.select或者标准ufunc的组合来表达实在不行再用frompyfunc或vectorize但要有心理准备它的速度不会比Python循环快太多只是代码简洁一些。顺带一提自定义ufunc返回的对象类型通常是object需要数值运算时记得转一下astype。4. 线性代数与随机数科学计算的左右手4.1 np.linalg模块速览先记住这几个函数NumPy的np.linalg模块提供了常见线性代数运算。先列个速查表函数作用典型场景np.linalg.inv矩阵求逆解线性系统、统计计算np.linalg.det行列式判断矩阵是否奇异、体积计算np.linalg.solve解线性方程组比先求逆再乘更稳定且更快np.linalg.eig特征值和特征向量PCA主成分分析np.linalg.svd奇异值分解降维、推荐系统、图像压缩np.linalg.norm向量/矩阵范数距离计算、正则化np.linalg.qrQR分解最小二乘、数值稳定的求特征值这里最值得强调的一点是解线性方程组Ax b时优先用np.linalg.solve(A, b)而不是np.linalg.inv(A) b。原因倒不是慢而是inv先显式计算了逆矩阵再乘b中间会引入额外的浮点误差而且计算量更大。solve内部用的是LU分解数值稳定性更好速度也更快。4.2 行列式与矩阵求逆的数值稳定性一个容易被忽略的坑np.linalg.det用来求行列式。行列式接近0意味着矩阵接近奇异也就是矩阵的某些行或列近乎线性相关这种情况下求逆的结果会在数值上非常不稳定。实际业务中行列式很难严格等于0往往是一个很小的值这时候inv的结果会出现巨大的数值波动。我之前处理过一个最小二乘问题设计矩阵的条件数特别大直接求逆得到的结果和用np.linalg.lstsq得到的结果差了三个数量级。查资料才发现线性最小二乘的“教科书公式”X^T X的逆再乘X^T在数值上是不稳定的。正确的做法是优先用np.linalg.lstsq或者np.linalg.pinv伪逆它们在内部会做SVD分解对接近奇异的情况处理得更好。所以给大家一个实操建议凡是涉及“求逆”的场景先想一想能不能换成solve、lstsq、pinv中的某一个。大多数情况下求逆只是数学推导里的中间步骤落到代码上并不需要真的算逆矩阵。4.3 np.random的进化从RandomState到Generator随机数生成也是NumPy中日常频率很高的部分。老版本的写法是np.random.seed(42)加np.random.randn(...)、np.random.randint(...)这套基于RandomState的接口虽然短但存在一些问题全局种子容易干扰其他模块的随机性而且随机数算法的选择不够灵活。NumPy 1.17开始推荐的新接口是Generator用法如下rng np.random.default_rng(42) # 自己的随机数发生器 x rng.normal(0, 1, size(3, 3)) # 标准正态分布 y rng.integers(0, 100, size10) # 整数随机数 idx rng.choice(100, size5, replaceFalse) # 不放回抽样Generator相比RandomState的几个好处随机数生成速度更快支持PCG64等更新更强的算法用起来更符合“每个过程一个随机源”的工程实践避免全局状态污染。如果你管理多个数据管线建议给每个阶段都初始化一个独立的default_rng而不是共用一个种子。随机数在实际中还能做很多有趣的事。比如用蒙特卡洛方法估算圆周率在一个边长为2的正方形内随机撒点统计落在内切圆内的比例乘以4就是圆周率的近似值rng np.random.default_rng(42) n 100_000 x rng.uniform(-1, 1, n) y rng.uniform(-1, 1, n) pi_approx 4 * ((x**2 y**2) 1).mean() print(pi_approx) # 约 3.14这段代码充分体现了“向量化随机数”组合的威力一个for循环都没有。5. 结构化数组与文件IO让ndarray带上语义5.1 结构化数组给每列起名字的ndarrayNumPy的ndarray默认是同构的也就是每个元素类型一致。但有一种情况例外结构化数组。它允许一个数组里包含多个字段每个字段有自己的名字和类型。可以把它理解成“更简洁的DataFrame雏形”。定义结构化数组非常简单只要在dtype里指定字段名和类型data np.array( [(1, 2.5, bapple), (2, 3.6, bbanana)], dtype[(id, i4), (score, f8), (name, S10)] ) print(data[score]) # [2.5 3.6] print(data[0][name]) # bapple有了字段名之后取整列数据就变成了data[score]不再需要记住“第2列是score”。这种代码在可读性上的提升非常大尤其是当列数多、字段语义重要的时候。结构化数组也可以和普通ndarray互相转换。如果你有一堆csv数据每列的语义各不相同用结构化数组存储会比用一个大的object数组清晰得多。不过要注意结构化数组的算力支持不如普通ndarray全面——很多ufunc是逐元素作用于整个结构体而不是按字段分别作用所以遇到复杂运算时通常还是建议把结构化数组拆成多个普通数组算完再合回结构化数组。5.2 np.loadtxt与np.genfromtxt文本数据读写的细节参数读写文本文件最常见的是np.loadtxt和np.savetxt。np.savetxt保存二维数组np.loadtxt读回来中间用fmt控制格式、delimiter控制分隔符arr np.array([[1.5, 2.5], [3.5, 4.5]]) np.savetxt(data.csv, arr, delimiter,, fmt%.6f) arr2 np.loadtxt(data.csv, delimiter,)loadtxt还有几个参数容易被忽略但很实用usecols只读某些列skiprows跳过文件头unpackTrue把结果按列拆开相当于转置输出。如果数据文件带有列名表头想直接按名字取列更合适的工具是np.genfromtxtarr np.genfromtxt(data_with_header.csv, delimiter,, namesTrue) print(arr[temperature])namesTrue表示第一行是列名读出来的结构化数组可以直接用字段名访问。genfromtxt还有一个好处支持处理缺失值missing_values和filling_values参数可以在读取阶段就把无效数据替换成指定值。5.3 np.savez与np.memmap处理超大数据集的两把钥匙当中间结果太多想一次性存下多个数组用np.savez或np.savez_compressed比逐个np.save方便得多np.savez_compressed(result.npz, XX, yy, labelslabels) loaded np.load(result.npz) X2 loaded[X].npz文件本质是一个zip包savez_compressed会做压缩对重复值多的数据效果尤其明显。再往外走一步如果数据大到根本装不进内存np.memmap就派上用场了。memmap把磁盘上的二进制文件映射到内存地址空间读写时操作系统按页加载不需要一次性读入整个数组。这在处理超大矩阵时非常关键# 创建一个 10000x10000 的 float64 文件约 800MB mm np.memmap(large.dat, dtypefloat64, modew, shape(10000, 10000)) # 可以像普通数组一样按块操作 mm[:, 0] 1 mm.flush()memmap的坑在于它默认返回的是一个memory-mapped数组很多操作会触发“物化”把整块数据读入内存比如np.array(mm)就会真正分配内存。所以用它时要有意识地做分块计算避免无意中把整个数组load进RAM。6. 视图vs副本与性能优化NumPy代码提速的几个实操点6.1 视图还是副本一个切片引发的“血案”对ndarray做基本切片arr[1:5]时返回的是原数组的视图。视图和数据是共享内存的所以修改视图会同步改动原数组。这对性能是好事避免了不必要的复制但对不熟悉的人来说这是一个经典的“血案”来源。a np.arange(10) b a[2:8] b[0] 999 print(a) # [ 0 1 999 3 4 5 6 7 8 9]上面这段代码里我只改了b结果a也跟着变了。很多新手写代码时以为切片是“复制了一份”结果无意间改了原数据排查半天。那么什么操作返回视图什么操作返回副本呢操作返回类型说明基础切片arr[1:5]视图与原数组共享内存整数/布尔索引arr[[0,2]]、arr[mask]副本新分配内存arr.reshape()视图可能在原数据连续时可共享内存arr.ravel()视图优先尽量返回视图arr.flatten()副本总是复制arr.copy()副本显式复制判断两个数组是否共享内存可以用np.shares_memory(a, b)。如果你确实需要一份独立的副本请显式调用.copy()这是最稳妥的做法。6.2 C-order与F-order为什么遍历顺序会影响性能NumPy数组在内存中的排列顺序有行优先C-order和列优先F-order之分。默认是C-order即先沿着行方向填满再换下一行F-order则先填满列方向。这个顺序对性能的影响在于CPU缓存和分页的局部性。如果你的算法是逐行遍历数组C-order会非常顺畅因为相邻元素在内存中也是相邻的如果算法是逐列遍历F-order会更快。用np.ascontiguousarray(arr)可以强制转换成C-ordernp.asfortranarray(arr)转成F-order。举一个典型的例子在按列累加时C-order的数组访问跨越大步长缓存命中率下降性能可能差好几倍。之前我做矩阵乘法优化时把其中一个矩阵转成F-order后耗时下降了约30%。所以如果你的算法有明确的遍历方向检查一下数组的arr.flags[C_CONTIGUOUS]或arr.flags[F_CONTIGUOUS]是值得养成习惯的一件事。6.3 实测提速案例用广播重构三层for循环最后用一个真实案例来收尾。假设我要计算两个点集A和B之间的平方距离矩阵A形状(500, 3)B形状(800, 3)。第一版用Python三层循环写严格说是两层循环因为每对点内部还要遍历3个坐标D np.empty((500, 800)) for i in range(500): for j in range(800): s 0 for k in range(3): s (A[i, k] - B[j, k]) ** 2 D[i, j] s第二版用广播加ufuncdiff A[:, None, :] - B[None, :, :] # (500, 800, 3) D (diff ** 2).sum(axis-1)实测下来前者在我的机器上大约耗时0.8秒左右后者只有约3毫秒速度提升了两百多倍。这个差距不是“稍微优化”的程度而是“明明可以秒开却要卡很久”的差别。还有一个中间的写法用scipy.spatial.distance.cdist它在底层调了向量化的C实现速度也很快代码更简洁。如果项目里已经装了scipy优先用专业工具也是完全合理的选择。不过NumPy版本的优点是不依赖额外库逻辑也足够直观适合作为理解距离计算的模板。最后再说一个我自己的习惯。写NumPy代码时我会先问自己三个问题这个操作能不能用广播写能不能用ufunc写是不是在无意中创建了副本带着这三个问题去看代码很多“慢”都能在没开始跑之前就规避掉。另外如果遇到UserWarning: failed to initialize numpy这类报错多半是环境里存在多个NumPy版本冲突或者某个第三方库引用了不同路径的NumPy。处理方式也很简单在干净环境里重装一遍NumPy或者用pip install --upgrade numpy统一版本。别慌这类问题大多不是代码逻辑的问题而是环境管理的问题。希望这篇笔记能帮你在NumPy上少走点弯路把省下来的时间留给真正需要思考的地方。