粒子物理学中不变质量的数值稳定计算？

Question

在粒子物理学中，我们必须大量计算invariant mass，这是用于二体衰变

当动量 (p1, p2) 与质量 (m1, m2) 相比有时非常大（高达 1000 倍或更多）时。在这种情况下，当在计算机上使用浮点数进行计算时，最后两项之间会发生很大的抵消。

可以使用什么样的数值技巧来准确计算任何输入？

问题是关于使用合适的数值技巧来提高浮点数计算的准确性，因此解决方案应该与语言无关。出于演示目的，首选 Python 中的实现。重新表述问题并增加基本运算量的解决方案是可以接受的，但建议使用其他数字类型（如十进制或多精度浮点数）的解决方案则不可接受。

注意：原始问题以 Python 表达式的形式提出了一个简化的 1D 维度问题，但问题是针对以 3D 维度给出动量的一般情况。问题以这种方式重新表述。

Answer 1

通过 Stackoverflow 上列出的一些技巧以及 Jakob Stark 在他的回答中描述的转换，可以将等式重写为不再遭受灾难性取消的形式。

原来的问题是求一维解，有一个简单的解，但在实践中，我们需要 3D 的公式，然后解就更复杂了。见this notebook for a full derivation。

Python中3D数值稳定计算的示例实现：

import numpy as np

# numerically stable implementation
@np.vectorize
def msq2(px1, py1, pz1, px2, py2, pz2, m1, m2):
    p1_sq = px1 ** 2 + py1 ** 2 + pz1 ** 2
    p2_sq = px2 ** 2 + py2 ** 2 + pz2 ** 2
    m1_sq = m1 ** 2
    m2_sq = m2 ** 2
    x1 = m1_sq / p1_sq
    x2 = m2_sq / p2_sq
    x = x1 + x2 + x1 * x2
    a = angle(px1, py1, pz1, px2, py2, pz2)
    cos_a = np.cos(a)
    if cos_a >= 0:
        y1 = (x + np.sin(a) ** 2) / (np.sqrt(x + 1) + cos_a) 
    else:
        y1 = -cos_a + np.sqrt(x + 1) 
    y2 = 2 * np.sqrt(p1_sq * p2_sq)
    return m1_sq + m2_sq + y1 * y2

# numerically stable calculation of angle
def angle(x1, y1, z1, x2, y2, z2):
    # cross product
    cx = y1 * z2 - y2 * z1
    cy = x1 * z2 - x2 * z1
    cz = x1 * y2 - x2 * y1
    
    # norm of cross product
    c = np.sqrt(cx * cx + cy * cy + cz * cz)
    
    # dot product
    d = x1 * x2 + y1 * y2 + z1 * z2
    
    return np.arctan2(c, d)

数值稳定的实现永远不会产生负面结果，这是幼稚实现的常见问题，即使是双精度也是如此。

让我们将数值稳定函数与一个简单的实现进行比较。

# naive implementation
def msq1(px1, py1, pz1, px2, py2, pz2, m1, m2):
    p1_sq = px1 ** 2 + py1 ** 2 + pz1 ** 2
    p2_sq = px2 ** 2 + py2 ** 2 + pz2 ** 2
    m1_sq = m1 ** 2
    m2_sq = m2 ** 2
    
    # energies of particles 1 and 2
    e1 = np.sqrt(p1_sq + m1_sq)
    e2 = np.sqrt(p2_sq + m2_sq)

    # dangerous cancelation in third term
    return m1_sq + m2_sq + 2 * (e1 * e2 - (px1 * px2 + py1 * py2 + pz1 * pz2))

对于下图，动量 p1 和 p2 是从 1 到 1e5 随机选取的，m1 和 m2 的值是从 1e-5 到 1e5 随机选取的。所有实现都以单精度获取输入值。两种情况下的引用都是使用 mpmath 使用带有 100 位小数的朴素公式计算得出的。

简单的实现会丢失某些输入的所有准确性，而数值稳定的实现不会。

Answer 2

如果你把例如表达式中的m1 = 1e-4、m2 = 1e-4、p1 = 1 和p2 = 1，你得到的4e-8 是双精度的，而0.0 是单精度计算的。我假设，您的问题是关于如何通过单精度计算获得4e-8。

您可以做的是上述表达式的泰勒展开式（大约 m1 = 0 和 m2 = 0）。

e ~ e|(m1=0,m2=0) + de/dm1|(m1=0,m2=0) * m1 + de/dm2|(m1=0,m2=0) * m2 + ...

如果我计算正确，零阶项和一阶项为 0，二阶展开为

e ~ (p1+p2)/p1 * m1**2 + (p1+p2)/p2 * m2**2

即使使用单精度计算，这也会准确地产生4e-8。如果需要，您当然可以在展开式中做更多项，直到达到单个浮点数的精度限制。

编辑
如果mi 并不总是比pi 小很多，你可以进一步按摩方程得到

复杂的部分现在是方括号中的部分。对于各种x 值，它本质上是sqrt(x+1)-1。如果x 非常小，我们可以使用平方根的泰勒展开式（例如here）。如果x 的值较大，则公式工作得很好，因为1 的加法和减法由于浮点精度而不再丢弃x 的值。所以x 的一个阈值必须选择低于一个切换到泰勒展开式。