Numpy float64矩阵行乘法出现数值不稳定问题如何解决？

问题原因 #

• 本质是浮点精度误差和底层BLAS库的运算调度共同导致：1/3是无限循环二进制小数，无法被float64类型精确存储，本身就存在微小的舍入误差。
• numpy的matmul底层调用优化后的BLAS接口，多线程场景下不同列的点积累加顺序不同，微小的误差在不同运算路径下就会出现你观察到的差异：前4列的累加残留了1e-17量级的误差，最后一列的累加顺序刚好抵消了该误差。
• 替换为1/2后结果稳定是因为1/2是2的整数次幂，可被float64精确表示，不存在舍入误差，无论运算顺序如何结果都一致。

解决办法 #

方案1：利用矩阵结构优化（最推荐，性能无损失） #

你的场景中x矩阵所有列取值完全相同，无需重复计算5次点积，仅计算1次后复制为对应长度的数组即可，完全避免不同运算路径带来的误差：

# 只计算第一列的点积，再复制为和x列数等长的数组
result = np.full(x.shape[1], np.dot(w, x[:, 0]))

方案2：误差截断 #

如果无法提前确定矩阵列重复，可根据业务精度要求设置合理的误差阈值，将小于阈值的微小误差统一置0：

# 阈值可根据实际场景调整，float64场景下通常设置为1e-15即可覆盖绝大多数舍入误差
error_threshold = 1e-15
result[np.abs(result) < error_threshold] = 0

方案3：强制单线程运算保证结果一致性 #

如果需要严格保证相同输入得到相同输出，可关闭BLAS的多线程调度，固定运算顺序，在导入numpy前设置环境变量即可：

import os
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np

注意该方案会降低大矩阵运算的性能。

方案4：使用更高精度的数值类型 #

如果对精度要求极高，可使用np.longdouble类型替换默认的np.float64，进一步降低舍入误差的影响：

x = np.zeros((4, 5,), dtype=np.longdouble)
w = np.array([1, 1, -1, 1 / 3], dtype=np.longdouble)

该方案的支持程度取决于运行平台的硬件和系统实现，性能也会略低于float64运算。

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