计算科学 - 如何可靠地添加大指数项而不会出现溢出错误？ - 吾爱随笔录

如何可靠地添加大指数项而不会出现溢出错误？

计算科学蒙特卡洛浮点统计数据

2021-12-06 21:00:58

马尔可夫链蒙特卡罗中一个非常常见的问题涉及计算大指数项之和的概率，

$e^{a_1} + e^{a_2} + ...$

其中的组件 $a$ 范围可以从非常小到非常大。我的方法是分解出最大的指数项 $K := \max_{i}(a_{i})$ 以便：

a^{'} = K + l o g (e^{a_{1} - K} + e^{a_{2} - K} + . . .)

$a' =K + log\left( e^{a_1 - K} + e^{a_2 - K } + ... \right)$

e^{a^{'}} \equiv e^{a_{1}} + e^{a_{2}} + . . .

$e^{a'} \equiv e^{a_1} + e^{a_2} + ...$

的所有元素 $a$ 都很大，则此方法是合理的，但如果不是，则不是一个好主意。当然，较小的元素无论如何都不会对浮点和做出贡献，但我不确定如何可靠地处理它们。在 R 代码中，我的方法如下所示：

if ( max(abs(a)) > max(a) )
  K <-  min(a)
else
  K <- max(a)
ans <- log(sum(exp(a-K))) + K

这似乎是一个足够普遍的问题，应该有一个标准的解决方案，但我不确定它是什么。感谢您的任何建议。

4个回答

有一个简单的解决方案，只需两次通过数据：

首先计算

K := max_{i} a_{i},

$K := \max_i\; a_i,$

它告诉你，如果有 $n$ 项，那么

\sum_{i} e^{a_{i}} \leq n e^{K} .

$\sum_i e^{a_i} \le n e^K.$

由于您可能没有 $n$ 接近甚至 $10^{20}$ 大，因此您不必担心在双精度计算

τ := \sum_{i} e^{a_{i} - K} \leq n

$\tau := \sum_i e^{a_i-K} \le n$ 时溢出.

因此，计算 $\tau$ 然后你的解决方案是 $e^K \tau$ 。

为了在将双打加在一起时保持精度，您需要使用Kahan Summation，这是相当于有一个进位寄存器的软件。

这对于大多数值来说都很好，但是如果您遇到溢出，那么您将达到IEEE 754 双精度的限制，大约是。此时您需要一个新的表示。您可以在加法时检测溢出 by ，也可以检测到大的指数以评估 by 。此时，您可以通过移动指数并跟踪此移动来修改双精度数的解释。 $e^{709.783}$ doubleMax - sumSoFar < valueToAddexponent > 709.783

这在很大程度上类似于您的指数偏移方法，但此版本以 2 为底，不需要初始搜索即可找到最大指数。因此。 $\mathit{value}\times2^\mathit{shift}$

#!/usr/bin/env python
from math import exp, log, ceil

doubleMAX = (1.0 + (1.0 - (2 ** -52))) * (2 ** (2 ** 10 - 1))

def KahanSumExp(expvalues):
  expvalues.sort() # gives precision improvement in certain cases 
  shift = 0 
  esum = 0.0 
  carry = 0.0 
  for exponent in expvalues:
    if exponent - shift * log(2) > 709.783:
      n = ceil((exponent - shift * log(2) - 709.783)/log(2))
      shift += n
      carry /= 2*n
      esum /= 2*n
    elif exponent - shift * log(2) < -708.396:
      n = floor((exponent - shift * log(2) - -708.396)/log(2))
      shift += n
      carry *= 2*n
      esum *= 2*n
    exponent -= shift * log(2)
    value = exp(exponent) - carry 
    if doubleMAX - esum < value:
      shift += 1
      esum /= 2
      value /= 2
    tmp = esum + value 
    carry = (tmp - esum) - value 
    esum = tmp
  return esum, shift

values = [10, 37, 34, 0.1, 0.0004, 34, 37.1, 37.2, 36.9, 709, 710, 711]
value, shift = KahanSumExp(values)
print "{0} x 2^{1}".format(value, shift)

你的方法是可靠的。

您不需要确切地知道，只要足够好就可以避免溢出。因此，您可以在进行任何 MCMC 采样之前通过分析 $K$ $K$

有一个 R 包可以快速有效地实现“log-sum-exp 技巧”

http://www.inside-r.org/packages/cran/matrixStats/docs/logSumExp

logSumExp 函数接受一个数字向量 lX 并输出 log(sum(exp(lX)))，同时使用您描述的方法避免下溢和上溢问题。

其它你可能感兴趣的问题

上一篇科学 Python 中的有限差分法推荐下一篇是否有找到渐近斜率的数值算法？