如何处理科学计算中的下溢？

Question

我正在处理概率模型，当对这些模型进行推断时，估计的概率可能会变得非常小。为了避免下溢，我目前正在log域中工作（将概率的log值存储）。概率的乘法等价于加法，并且使用以下公式进行求和：

log(exp(a) + exp(b)) = log(exp(a - m) + exp(b - m)) + m

假设 m = max(a, b)。

我使用了一些非常大的矩阵，并且必须对这些矩阵进行逐元素指数运算以计算矩阵-向量乘积。这一步骤非常费时，我想知道在处理概率时是否存在其他方法来处理下溢。

编辑：出于效率考虑，我正在寻找一种使用原始类型而不是存储实数任意精度表示的对象的解决方案。

编辑2：我正在寻找比log域技巧更快的解决方案，而不是更准确的解决方案。我对当前获得的精度感到满意，但我需要更快的方法。特别是，在矩阵-向量乘积期间会发生求和，并且我希望能够使用高效的BLAS方法。

解决方案：在与Jonathan Dursi的讨论后，我决定将每个矩阵和向量因子分解为其最大元素，并将该因子存储在对数域中。乘法很简单。在加法之前，我必须通过两个因子的比率因子分解一个要添加的矩阵/向量。我每十个操作更新一次因子。

- Edouard

---

你必须使用Java吗？还是可以使用其他编程语言？ - enzom83

@Peter - 这并不是不寻常的情况。例如，在使用极大似然估计时，看到这样的数字并不罕见。即使起始点不尽如人意，您的优化器仍必须能够收敛。如果在那里出现下溢，则无法收敛。 - user85109

它似乎是一个相当抽象的问题。如果你用普朗克单位来衡量宇宙的年龄，你会得到大约2e58个时间单位，也就是任何事情可能发生的时间单位数量。如果某件事的概率小于1e-300，很难想象它不是相当接近不可能或者说理论上无法测量和不可知的。想一想你需要采取多少次测量才能知道某件事的概率是1e-58。 - Peter Lawrey

@Peter - 假设你正在对一条线上移动的粒子进行建模，其行为如下：在每个时间步长中，它可以向前移动一步，概率为0.5，或向后移动一步，概率也为0.5。长度为1000的一系列位置序列的概率为0.5^1000。通过一次测量，我观察到了一个概率非常低的序列。 - Edouard

只是一个提示：a-m或b-m中的一个将为0，因此您可以将1放在最大值的位置。 - kloop

Answer 1

This issue has come up recently on the 计算科学堆栈交换站 as well, and although there the immediate worry there was overflow, the issues are more or less the same.

Transforming into log space is certainly one reasonable approach. Whatever space you're in, to do a large number of sums correctly, there's a couple of methods you can use to improve the accuracy of your summations. Compensated summation approaches, most famously Kahan summation, keep both a sum and what's effectively a "remainder"; it gives you some of the advantages of using higher precision arithmetic without all of the cost (and only using primitive types). The remainder term also gives you some indication of how well you're doing.

除了改进你的加法机制之外，改变添加项的顺序可以产生很大的差异。将项排序，使得从最小到最大的求和可以帮助，因为这样你不再频繁地添加非常不同的项（这可能会导致重大的舍入问题）；在某些情况下，进行log₂N次重复的成对求和也比直接进行线性求和更好，具体取决于你的项的形式。

所有这些方法的有用性都在很大程度上取决于数据的属性。虽然任意精度数学库使用起来计算时间（和可能的内存）非常昂贵，但它们具有相当通用的解决方案的优点。

- Jonathan Dursi

---

谢谢你提供非常有趣的答案。然而，我正在寻找一种更有效的方法，而不是更准确的方法（我对使用对数域技巧得到的准确性感到满意）。而且，如果不在对数空间中工作，仅使用补偿求和只能解决准确性问题，而不能解决下溢风险。 - Edouard

你不关心准确性，但是担心下溢吗？下溢不是准确性的考虑因素吗？我不认为我理解你在寻找什么。 - Jonathan Dursi

我所说的“准确性”是指求和的准确性。使用补偿求和，即使在乘以两个小数时也可以获得无法用“double”表示的数字。在对长HMM进行推断时，您可以获得比“10 ^ -324”更小但具有相同数量级的中间量。通过最大值进行因式分解，可以计算出准确的总和。这就是我的当前解决方案。基本上，我正在寻找一种表示小数的方法，具有高效的加法和乘法。现在我只有高效的乘法。 - Edouard

期望的指数范围和结果的数字精度是多少？ - Jonathan Dursi

我会尝试将这个和我的旧解决方案结合起来。每个向量或矩阵都将通过其最大元素进行因式分解，并存储在对数域中。乘法和加法很简单。我只需要每十到二十次操作更新一次因子项。感谢您的帮助。 - Edouard

Answer 2

选项1：Commons Math - The Apache Commons Mathematics Library

Commons Math是一个轻量级的自包含数学和统计组件库，用于解决Java编程语言或Commons Lang中最常见的问题。

注意：API保护构造函数以强制使用工厂模式，同时将工厂命名为DfpField（而不是更直观的DfpFac或DfpFactory）。因此，您必须使用

new DfpField(numberOfDigits).newDfp(myNormalNumber)

要实例化一个Dfp，然后可以在此上调用.multiply或其他内容。我想提一下这一点，因为它有点令人困惑。

选项2：GNU Scientific Library或Boost C++ Libraries。在这些情况下，您应该使用JNI来调用这些本地库。

选项3：如果您可以自由使用其他程序和/或语言，则可以考虑使用用于数值计算的程序/语言，例如Octave、Scilab等。

选项4：Java的BigDecimal。

- enzom83

---

至少Matlab和Octave都有一些Java绑定。 - Voo

谢谢提供的参考资料，但我认为它们对我没用。选项1和4：使用任意精度十进制数太昂贵了，因为它们使用对象而不是原始类型，并且使用这种表示法进行加法和乘法计算更加昂贵。选项2：与1和4相同的问题（AFAIK），并且我更喜欢坚持使用Java。选项3：我已经使用NumPy和Matlab有一段时间了，出现了相同的问题，因为它们也使用浮点数和双精度浮点数。 - Edouard

@Edouard：但从这个角度来看，Java是最不适合模拟的语言，因为它是一种“半编译”（即“半解释”）语言，所以你会遇到性能问题。相反，Octave、Scilab和类似的语言都有自己优化矩阵和向量操作的例程，事实上它们经常用于模拟。然而我记得在Matlab中你可以设置精度：看看这个链接。 - enzom83

@enzom83 - 我已经使用这些工具三年了（主要是 Scilab 和 Numpy）。当我在小型隐马尔可夫模型上进行推理时，我使用了广为人知的对数域技巧。但即便是对于这些语言，指数计算仍然是瓶颈。 - Edouard

@enzom83 - 我正在使用与BLAS绑定来进行矩阵计算，因此矩阵/向量操作的性能在这里不是问题。而且相比C++甚至FORTRAN，Java也不差：请参见language shootout。 - Edouard

Answer 3

我几年前遇到了类似的问题。解决方法是开发出一个log（1 + exp（-x））的近似值。该近似值的范围并不需要很大（x从0到40就足够了），至少在我这种情况下，精度也不需要特别高。

在您的情况下，看起来您需要计算log（1 + exp（-x1）+ exp（-x2）+ ...）。舍弃那些很大的负值。例如，假设a，b和c是三个对数概率，并且0> a> b> c。如果a-c> 38，则可以忽略c。它根本不会对您的联合对数概率产生贡献，至少如果您正在使用双精度时不会有贡献。

- David Hammen

---

聪明的技巧。但我认为开发一个比对n个双精度数取exp函数更快的log(1 + exp(x1) + exp(x2) + ...)的近似值是相当具有挑战性的。 - Edouard

你仍然可以使用排除那些极低概率事件的技巧。如果你正在使用IEEE双精度浮点数，1+exp(-37)恰好等于1。这将立即解决你的下溢问题。 - David Hammen

Answer 4

与其以对数形式存储值，我认为您最好使用与double相同的概念，即浮点表示。例如，您可以将每个值存储为两个long，一个用于符号和尾数，另一个用于指数。（真实浮点具有经过精心调整的设计，以支持许多边缘情况并避免浪费单个位；但您可能不需要太担心这些问题，并且可以专注于以简单易行的方式设计它。）

- ruakh

---

OP正在研究概率模型。在这类问题中，对数概率非常常见。 - David Hammen

我考虑过这个问题。但是正如我在编辑后的问题中所说，出于效率原因，我更喜欢坚持使用基本类型（doubles），而不是开发一种更适合我的需求但会导致性能问题的新类型。 - Edouard

@Edouard：我不知道。对我来说，使用两个long并执行普通的整数算术比使用一个double并执行对数和指数运算更慢似乎很奇怪，但我会相信你的话。 - ruakh

我表达不够清楚。我的意思是，使用自定义类型会迫使我编写线性代数函数来处理该类型的矩阵/向量计算。但是与 BLAS 或其他针对“double”进行优化的线性代数包竞争有点困难。 - Edouard

Answer 5

我不明白为什么这个公式可行，但是这个公式似乎可行且更简单：

c = a + log(1 + exp(b - a))

其中 c = log(exp(a)+exp(b))。

- Houshalter