计算sqrt(x+a) - sqrt(x)时如何保证数值稳定性？

Question

计算sqrt(x+a) - sqrt(x)时如何保证数值稳定性？

9

在x和a >= 0的全部参数范围内，是否有一种优雅的方法来进行数值稳定的评估以下表达式？

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

是否有任何编程语言或库提供此类功能？如果有，它的名称是什么？我现在没有使用上面的表达式的具体问题，但过去遇到过很多次，并且一直认为这个问题肯定已经解决了！

- Andreas H.

3

一些库，尤其是Boost，提供了一个名为sqrt1pm1()的函数，旨在准确计算sqrt(x+1)-1。如果您已经使用了这样的库，您可以使用该函数以数值稳定的方式实现sqrt(x+a)-sqrt(x)，即sqrt1pm1(a/x)*sqrt(x)。请注意，在实现中不要改变原来的意思。 - njuffa

@njuffa：啊，非常有趣。虽然log1p和expm1这样的函数很常见，但我从未遇到过sqrt1pm1。一方面，当我们可以很容易地仿真这个功能时，单独创建一个函数似乎有些奇怪。另一方面，如果C标准库中有这个函数，我肯定会经常用到它。 - Mark Dickinson

@MarkDickinson 正如Kahan所示，log1p和expm1也很容易模拟。假定在库中提供这样的函数的目的是为了向不太了解数值分析的程序员提供最快速和最准确的实现。 - njuffa

@njuffa：我认为这是一个非常不同的“容易”的价值观。例如，请参阅Python源代码，了解在尝试模拟log1p时涉及到的考虑因素，特别是在既不能依赖浮点格式也不能依赖当前舍入模式的情况下。 - Mark Dickinson

@njuffa：了解sqrtpm1函数的信息很有趣。这实际上是我想学习的内容...该方法的缺点当然是必须添加x=0的特殊情况... - Andreas H.

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Mark Dickinson · Accepted Answer

是的，可以！只要 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。让我们尝试一个特定的小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双精度浮点数的上下文中，立方和更高阶项是微不足道的，该标准只提供大约16个十进制位的精度。为了比较，让我们计算一下该值：

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