Python中np.divide(x,y)与x/y结果差异的原因探究

最近在代码中发现一个问题：将np.divide(x, y)替换为x / y后，之前的异常结果得以修复。此前根据NumPy文档描述，认为np.divide(x, y)与x / y等价，但实际运行结果存在明显差异。问题已解决，现说明背后原因：

测试代码 #

import numpy as np


x1 = np.array([[281], [15831], [30280], [975], [313], [739], [252], [10364], [21480], [1447], [315], [772], [95], [2710], [7408], [215], [111], [158], [0], [88], [21], [661], [0], [0], [0], [5], [4], [0], [12], [0], [0], [50], [28], [0], [0], [272]])
x2 = np.array([[499], [6315], [33800], [580], [208], [464], [384], [3127], [19596], [2319], [218], [1740], [217], [411], [4250], [223], [406], [267], [2], [0], [16], [0], [0], [0], [0], [8], [3], [0], [18], [0], [1], [0], [41], [0], [0], [0]])
x3 = np.array([[507], [6180], [34005], [555], [200], [451], [390], [3024], [19492], [2425], [211], [1848], [223], [396], [4097], [224], [406], [282], [2], [0], [16], [0], [0], [0], [0], [8], [3], [0], [19], [0], [2], [0], [45], [0], [0], [0]])
x4 = np.array([[507], [6178], [34017], [554], [200], [451], [391], [3022], [19486], [2439], [210], [1865], [223], [396], [4089], [224], [406], [284], [2], [0], [16], [0], [0], [0], [0], [8], [3], [0], [19], [0], [2], [0], [46], [0], [0], [0]])

not_zero = (x1 + x2) != 0
x = np.divide(2*(x1 - x2)**2, x1 + x2, where=not_zero)
r = (2*(x1[not_zero] - x2[not_zero])**2) / (x1[not_zero] + x2[not_zero])
print("n1 =",x.max(),"	t1 =", r.max())

not_zero = (x2 + x3) != 0
x = np.divide(2*(x2 - x3)**2, x2 + x3, where=not_zero)
r = (2*(x2[not_zero] - x3[not_zero])**2) / (x2[not_zero] + x3[not_zero])
print("n2 =",x.max(),"	t2 =", r.max())

not_zero = (x3 + x4) != 0
x = np.divide(2*(x3 - x4)**2, x3 + x4, where=not_zero)
r = (2*(x3[not_zero] - x4[not_zero])**2) / (x3[not_zero] + x4[not_zero])
print("n3 =",x.max(),"	t3 =", r.max())

运行输出 #

n1 = 8177.933351395286  t1 = 8177.933351395286
n2 = 873842.0           t2 = 6.501672240802676
n3 = 1322.0             t3 = 0.15566927013196877

运行环境 #

• Python版本：3.7.6
• NumPy版本：1.17.0

原因分析 #

这是NumPy 1.17.0版本中的已知bug，在后续更新版本（如1.18.0及以上）中已修复。

早期NumPy版本中，ufunc（通用函数，如np.divide）的where参数实现存在缺陷：当使用where指定计算条件时，部分满足条件的位置未正确执行除法操作，反而保留了被除数的原始值（例如2*(x2-x3)**2的数值），导致结果异常。

而直接使用x/y的方式，是先通过索引过滤出满足条件的元素再计算，完全避开了np.divide中where参数的实现问题，因此结果正确。若升级NumPy到修复后的版本，重新运行代码会发现np.divide与x/y的结果完全一致。

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