sqrt(x+a) - sqrt(x) 的数值稳定评估

Question

对于整个参数范围 x,a >= 0，是否有一种优雅的方法可以在数值上稳定地评估以下表达式？

f(x,a) = sqrt(x+a) - sqrt(x)

还有任何编程语言或库提供这种功能吗？如果有，用什么名字？我现在使用上面的表达式没有具体问题，但过去遇到过很多次，一直认为这个问题以前一定已经解决了！

Answer 1

是的，有！前提是x和a中至少有一个为正，可以用：

f(x, a) = a / (sqrt(x + a) + sqrt(x))

它在数值上是完全稳定的，但它本身几乎不值得一个库函数。当然，当x = a = 0时，结果应该是0.

解释：sqrt(x + a) - sqrt(x)等于(sqrt(x + a) - sqrt(x)) * (sqrt(x + a) + sqrt(x)) / (sqrt(x + a) + sqrt(x))。现在将前两项相乘得到 sqrt(x+a)^2 - sqrt(x)^2，它简化为 a.

这是一个证明稳定性的例子：原始表达式的麻烦情况是 x + a 和 x 的值非常接近（或者等效地，当 a 的值小得多时幅度比 x)。例如，如果 x = 1 和 a 很小，我们从围绕 1 的泰勒展开得知 sqrt(1 + a) 应该是 1 + a/2 - a^2/8 + O(a^3)，所以 sqrt(1 + a) - sqrt(1)应该接近 a/2 - a^2/8。让我们尝试选择 small a。这是原始函数（在本例中用 Python 编写，但您可以将其视为伪代码）：

def f(x, a):
    return sqrt(x + a) - sqrt(x)

这里是稳定版：

def g(x, a):
    if a == 0:
        return 0.0
    else:
        return a / ((sqrt(x + a) + sqrt(x))

现在让我们看看 x = 1 和 a = 2e-10 得到了什么：

>>> a = 2e-10
>>> f(1, a)
1.000000082740371e-10
>>> g(1, a)
9.999999999500001e-11

我们应该得到的值是（达到机器精度）：a/2 - a^2/8 - 对于这个特定的 a，三次和高阶项在 IEEE 754 double- 的上下文中是微不足道的精度浮点数，它只提供大约 16 位十进制数字的精度。让我们计算该值以进行比较：

>>> a/2 - a**2/8
9.999999999500001e-11

sqrt(x+a) - sqrt(x) 的数值稳定评估

Numerically stable evaluation of sqrt(x+a) - sqrt(x)

math

floating-point

numeric