处理含极小值的大型NumPy数组时的精度丢失问题

处理三个大型NumPy数组（元素为极小值）的计算时，单独计算数组元素结果符合预期，但对整个数组运算时出现精度丢失。相关变量particle_px、particle_py、particle_pz长度一致，particle_px与particle_py元素数值相近，particle_pz数值极小可忽略。

问题代码 #

def calculate_particle_energy(self):
    """
    calculates the particle energes from the momenta in the current file. output in keV
    """
    print("particle_px is of " + str(type(self.particle_px)) + " with " + str(len(self.particle_px)) + " elements")
    print("particle_px[0] = " + str(self.particle_px[0]))
    a = (self.particle_px**2 + self.particle_py**2 + self.particle_pz**2) / (2 * self.particle_mass*1.6022*(10**-16))
    print("denominator = " + str(2 * self.particle_mass*1.6022*(10**-16)))
    print("numerator_single = " + str(self.particle_px[0]**2 + self.particle_py[0]**2 + self.particle_pz[0]**2))
    print("numerator_array = " + str((self.particle_px**2 + self.particle_py**2 + self.particle_pz**2)[0]))
    print("correct a0 = " + str((self.particle_px[0]**2 + self.particle_py[0]**2 + self.particle_pz[0]**2) / (2 * self.particle_mass*1.6022*(10**-16))))
    print("wrong a0 = " + str(a[0]))

集群运行输出（NumPy 1.18.1） #

particle_px is of <class 'numpy.ndarray'> with 48933106 elements
particle_px[0] = 1.1980471e-25
denominator = 2.9190097272000005e-46
numerator_single = 1.817622024655414e-50
numerator_array = 0.0
correct a0 = 6.226844699140246e-05
wrong a0 = nan

本地测试（NumPy 1.24.1，小数组） #

测试代码 #

import numpy as np

class particle_energy:

    def __init__(self):
        self.particle_px = np.array([2e-26,4e-26,5e-26])
        self.particle_py = np.array([2e-26,4e-26,5e-26])
        self.particle_pz = np.array([2e-26,4e-26,5e-26])
        self.particle_mass = 9.10938 * (10 ** -31)
        
    def calculate_particle_energy(self):
        """
        calculates the particle energes from the momenta in the current file. output in keV
        """
        print("particle_px is of " + str(type(self.particle_px)) + " with " + str(len(self.particle_px)) + " elements")
        print("particle_px[0] = " + str(self.particle_px[0]))
        a = (self.particle_px**2 + self.particle_py**2 + self.particle_pz**2) / (2 * self.particle_mass*1.6022*(10**-16))
        print("denominator = " + str(2 * self.particle_mass*1.6022*(10**-16)))
        print("numerator_single = " + str(self.particle_px[0]**2 + self.particle_py[0]**2 + self.particle_pz[0]**2))
        print("numerator_array = " + str((self.particle_px**2 + self.particle_py**2 + self.particle_pz**2)[0]))
        print("correct a0 = " + str((self.particle_px[0]**2 + self.particle_py[0]**2 + self.particle_pz[0]**2) / (2 * self.particle_mass*1.6022*(10**-16))))
        print("wrong a0 = " + str(a[0]))

def main():
    pen = particle_energy()
    pen.calculate_particle_energy()
    
if __name__ == '__main__':
    main()   

测试输出 #

particle_px is of <class 'numpy.ndarray'> with 3 elements
particle_px[0] = 2e-26
denominator = 2.9190097272000005e-46
numerator_single = 1.2e-51
numerator_array = 1.2e-51
correct a0 = 4.1109832174183095e-06
wrong a0 = 4.1109832174183095e-06

提出的问题 #

1. 假设是处理大型数组时出现精度丢失，该假设是否正确？
2. 若正确，如何避免精度丢失并得到正确计算结果？
3. 兼顾精度与速度的最优计算方式是什么？

1. 精度丢失假设的正确性 #

不完全准确，核心原因并非数组大小，而是NumPy版本差异+数组存储 dtype 限制：

• 老版本NumPy（1.18.1）处理大规模数组时，若数组 dtype 为单精度浮点数（float32），极小值的平方会超出单精度最小可表示范围（约1.175e-38），直接下溢为0；但单独访问元素计算时，Python会自动将其提升为双精度（float64）运算，因此能得到正确结果。
• 本地新版本NumPy（1.24.1）对小数值运算的精度处理更优化，且测试用小数组默认使用float64 dtype，因此未复现问题。

2. 避免精度丢失的解决方法 #

• 转换数组 dtype：将数组从float32转为float64，双精度最小正浮点数约为2.225e-308，完全能容纳极小值的平方结果：

  self.particle_px = self.particle_px.astype(np.float64)
  self.particle_py = self.particle_py.astype(np.float64)
  self.particle_pz = self.particle_pz.astype(np.float64)
  

• 数值缩放法：给所有动量值乘以缩放因子（如1e25），将极小值转换为常规范围数值后计算，最后再缩放回原单位：

  scale_factor = 1e25
  px_scaled = self.particle_px * scale_factor
  py_scaled = self.particle_py * scale_factor
  pz_scaled = self.particle_pz * scale_factor
  numerator_scaled = px_scaled**2 + py_scaled**2 + pz_scaled**2
  numerator = numerator_scaled / (scale_factor**2)
  a = numerator / (2 * self.particle_mass*1.6022*(10**-16))
  

• 升级NumPy版本：若集群允许，升级到较新版本，其对浮点数运算的精度兼容性更好。

3. 兼顾精度与速度的最优方式 #

优先选择转换为float64 dtype：

• 精度上：双精度完全满足极小值运算需求，不会出现下溢；
• 速度上：现代CPU对float64运算支持成熟，大规模数组运算的性能损失极小，远低于其他精度补偿方案的额外开销。
  若集群环境限制无法转换 dtype，可采用数值缩放法，该方法不改变数组存储类型，仅通过算术缩放避免下溢，运算速度几乎与原数组一致。

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