np.hypot与np.subtract.outer为何比NumPy普通广播更快？有无更快的距离矩阵计算方法？

先给你拆解清楚这两个现象背后的原因，再推荐几个效率更高的实现方式：

一、为什么np.subtract.outer比普通广播快？ #

你测试里的普通广播写法obs[:, np.newaxis, :] - interp会先生成一个三维数组（形状(50000, 30000, 2)），这个数组光是存储就要占50000*30000*2*8 = 24GB的内存——这么大的数组很容易触发系统的虚拟内存交换（把数据写到硬盘），速度直接暴跌。

而np.subtract.outer是分维度单独计算，生成两个二维数组（d0和d1，每个形状(50000, 30000)），总内存其实和三维数组一样，但胜在：

1. 内存访问更高效：二维数组的内存布局更连续，CPU缓存命中率更高，避免了三维数组带来的缓存颠簸；
2. 底层优化更到位：subtract.outer本质是调用了BLAS（基础线性代数子程序）里的高效减法实现，而普通广播的三维减法没有用到这么针对性的优化，更多是通用的逐元素操作。

二、为什么np.hypot比sqrt(d0²+d1²)快？ #

手动计算sqrt(d0**2 + d1**2)会产生三个临时数组：d0²、d1²、两者的和，每一步都要读写大内存数组，而内存带宽恰恰是这种大数组计算的瓶颈。

反观np.hypot：

1. 无中间临时数组：它直接读取d0和d1的对应元素，原地计算sqrt(a² + b²)后直接写入结果数组，完全跳过了存储平方和、求和结果的步骤，大幅减少内存读写；
2. 底层指令优化：hypot是专门为计算斜边设计的函数，底层用了SSE/AVX等向量指令集加速，甚至会处理数值溢出的情况（比如超大数平方导致的精度丢失），比手动拼接的操作效率高得多。

三、更快的距离矩阵计算方法 #

如果你追求极致效率，这两个方法比你现有的实现快得多：

方法1：用scipy.spatial.distance.cdist #

这是专门为计算两组点距离矩阵设计的工具，底层是高度优化的C代码，还支持多线程（如果你的SciPy编译时开启了OpenMP），代码也极简：

from scipy.spatial.distance import cdist

# 直接计算欧氏距离矩阵
dist_matrix = cdist(obs, interp, metric='euclidean')

对于你的测试规模，这个方法的速度会比outer+hypot快2-5倍，甚至更多。

方法2：基于矩阵乘法的欧氏距离公式 #

欧氏距离的平方可以展开为：||x-y||² = ||x||² + ||y||² - 2x·y，利用这个公式可以用BLAS优化的矩阵乘法来计算，速度爆炸：

%%time
# 计算每个点的范数平方
obs_norms = np.sum(obs**2, axis=1)[:, np.newaxis]  # 形状(50000,1)
interp_norms = np.sum(interp**2, axis=1)[np.newaxis, :]  # 形状(1,30000)
# 计算点积矩阵（BLAS优化，速度极快）
dot_product = np.dot(obs, interp.T)
# 计算距离平方，再开根号
dist_sq = obs_norms + interp_norms - 2 * dot_product
dist_matrix = np.sqrt(dist_sq)

这个方法的核心是np.dot调用了BLAS里的gemm（通用矩阵乘法），这是CPU最擅长的优化操作，比逐元素的减法、平方快一个数量级，而且内存占用也更小（不需要存储d0和d1两个大数组）。

内容的提问来源于stack exchange，提问作者eroot163pi