在x和a >= 0的全部参数范围内,是否有一种优雅的方法来进行数值稳定的评估以下表达式?
f(x,a) = sqrt(x+a) - sqrt(x)
是否有任何编程语言或库提供此类功能?如果有,它的名称是什么?我现在没有使用上面的表达式的具体问题,但过去遇到过很多次,并且一直认为这个问题肯定已经解决了!
在x和a >= 0的全部参数范围内,是否有一种优雅的方法来进行数值稳定的评估以下表达式?
f(x,a) = sqrt(x+a) - sqrt(x)
是否有任何编程语言或库提供此类功能?如果有,它的名称是什么?我现在没有使用上面的表达式的具体问题,但过去遇到过很多次,并且一直认为这个问题肯定已经解决了!
是的,可以!只要 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
a=0
编辑了一个特殊情况。感谢您的纠正! - Mark Dickinsonx
和/或 a
是次正规数,同时继续假设 x
和 a
都是非负数,那么它们的平方根仍然在正常范围内,所以我看不出你会失去超过几个最不精确的单位。 (通过一些时间和足够大的纸张,应该可以得出实际上限误差的估值。)但是,现在我注意到我没有坚持 x
和 a
都是非负数。需要再思考一下。 - Mark Dickinson
sqrt1pm1()
的函数,旨在准确计算sqrt(x+1)-1。如果您已经使用了这样的库,您可以使用该函数以数值稳定的方式实现sqrt(x+a)-sqrt(x)
,即sqrt1pm1(a/x)*sqrt(x)
。请注意,在实现中不要改变原来的意思。 - njuffalog1p
和expm1
这样的函数很常见,但我从未遇到过sqrt1pm1
。一方面,当我们可以很容易地仿真这个功能时,单独创建一个函数似乎有些奇怪。另一方面,如果C标准库中有这个函数,我肯定会经常用到它。 - Mark Dickinsonlog1p
和expm1
也很容易模拟。假定在库中提供这样的函数的目的是为了向不太了解数值分析的程序员提供最快速和最准确的实现。 - njuffalog1p
时涉及到的考虑因素,特别是在既不能依赖浮点格式也不能依赖当前舍入模式的情况下。 - Mark Dickinson