6

如果我想取浮点数列表的乘积,那么通过添加它们的对数然后取总和的 exp 而不是仅仅将它们相乘,最坏情况/平均情况的精度是多少。有没有这种情况实际上更精确?

4

1 回答 1

10

没有任何上溢或下溢恶作剧,如果ab是浮点数,则乘积a*b将被计算在 1/2 ulp 的相对误差范围内。

因此,在乘以N doubles 链后,相对误差的粗略界限导致答案最多为 (1 - epsilon/2) -N的一个因子,约为 exp(epsilon N/2)。我想N在平均情况下,您可以预期 epsilon sqrt( ) 左右的偏差。(首先,这大约是 N epsilon。)

但是,这种策略更有可能发生指数溢出和下溢;由于次正规数的舍入,您更有可能得到无穷大、零和 NaN 以及不精确的值。

从这个意义上说,另一种方法更健壮,但在直接方法不会导致上溢或下溢的情况下,它会慢得多并且错误更严重。在 N 至少比 2 53小几个数量级的情况下,这是对标准双精度的非常非常粗略的分析:

你总是可以取一个有限浮点数的对数并得到一个有限浮点数,所以我们很酷。您可以直接将N浮点数相加以获得Nepsilon 最坏情况“相对”错误和 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 的差异。我敢肯定还有更糟糕的例子。

于 2013-04-07T03:45:23.140 回答