记录日志并添加与乘法

Question

如果我想获取浮点数列表的乘积,那么最坏情况/平均情况下的精度损失是通过添加它们的日志然后获取总和而不是仅仅乘以它们来实现的.有没有这种情况实际上更准确？

Answer 1

不存在任何溢或下溢的恶作剧,如果a和b是浮点数,则产品a*b将1/2 ULP的相对误差内来计算到.

N double因此,在乘以s 链之后对相对误差的粗略约束导致得到最多(1-epsilon/2)^-N的因子,其约为exp(ε N/ 2).我想你可以预期N平均情况下epsilon sqrt()的偏差.(首先,这是关于N epsilon.)

但是,该策略更有可能发生指数溢出和下溢; 由于次正规的四舍五入,你更有可能得到无穷大,零和NaN以及不精确的值.

另一种方法在这个意义上更加强大,但在直接方法不会导致溢出或下溢的情况下,它会慢得多,而且更糟糕.这是对标准双打的一个非常非常粗略的分析,其中N至少比2 ⁵³小几个数量级:

您总是可以获取有限浮点数的对数并获得有限的浮点数,因此我们在那里很酷.您可以直接添加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的差异.我敢肯定有更糟糕的例子.