如果我想取浮点数列表的乘积,那么通过添加它们的对数然后取总和的 exp 而不是仅仅将它们相乘,最坏情况/平均情况的精度是多少。有没有这种情况实际上更精确?
1 回答
没有任何上溢或下溢恶作剧,如果a
和b
是浮点数,则乘积a*b
将被计算在 1/2 ulp 的相对误差范围内。
因此,在乘以N
double
s 链后,相对误差的粗略界限导致答案最多为 (1 - epsilon/2) -N的一个因子,约为 exp(epsilon N
/2)。我想N
在平均情况下,您可以预期 epsilon sqrt( ) 左右的偏差。(首先,这大约是 N epsilon。)
但是,这种策略更有可能发生指数溢出和下溢;由于次正规数的舍入,您更有可能得到无穷大、零和 NaN 以及不精确的值。
从这个意义上说,另一种方法更健壮,但在直接方法不会导致上溢或下溢的情况下,它会慢得多并且错误更严重。在 N 至少比 2 53小几个数量级的情况下,这是对标准双精度的非常非常粗略的分析:
你总是可以取一个有限浮点数的对数并得到一个有限浮点数,所以我们很酷。您可以直接将N
浮点数相加以获得N
epsilon 最坏情况“相对”错误和 sqrt(N) epsilon 预期“相对”错误,或者使用Kahan 求和获得大约 3 epsilon 最坏情况“相对”错误。吓人的引号围绕“相对”,因为误差与您要求和的事物的绝对值的总和有关。
请注意,没有有限double
的对数的绝对值大于 710 左右。这意味着我们使用 Kahan 求和计算的对数和的绝对误差最多为 2130 N epsilon。当我们对对数和求幂时,我们从正确答案中得到最多 exp(2130 N epsilon) 的因子。
log-sum-exp 方法的病态示例:
int main() {
double foo[] = {0x1.000000000018cp1023, 0x1.0000000000072p-1023};
double prod = 1;
double sumlogs = 0;
for (int i = 0; i < sizeof(foo) / sizeof(*foo); i++) {
prod *= foo[i];
sumlogs += log(foo[i]);
}
printf("%a %a\n", foo[0], foo[1]);
printf("%a %a %a\n", prod, exp(sumlogs), prod - exp(sumlogs));
}
在我的平台上,我得到了 0x1.fep-44 的差异。我敢肯定还有更糟糕的例子。