对数求和技巧为什么不适合递归实现？

Question

我一直在研究log-sum-exp问题。我有一个以对数形式存储的数字列表，我想将它们相加并以对数形式存储。

朴素算法如下：

def naive(listOfLogs):
    return math.log10(sum(10**x for x in listOfLogs))

许多网站，包括： C语言中的logsumexp实现? 和 http://machineintelligence.tumblr.com/post/4998477107/ 推荐使用。

def recommend(listOfLogs):
    maxLog = max(listOfLogs)
    return maxLog + math.log10(sum(10**(x-maxLog) for x in listOfLogs))

即

def recommend(listOfLogs):
    maxLog = max(listOfLogs)
    return maxLog + naive((x-maxLog) for x in listOfLogs)

我不明白的是，如果推荐算法更好，为什么要递归调用它？这样做能带来更多的好处吗？

def recursive(listOfLogs):
    maxLog = max(listOfLogs)
    return maxLog + recursive((x-maxLog) for x in listOfLogs)

我的问题是：在这个计算中是否有其他技巧可以使其更加数值稳定？

- Jacob

---

我刚刚发现了scipy.misc.logsumexp：http://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.logsumexp.html - isti_spl

Answer 1

为了给其他人提供背景信息：当您直接计算以下类型的表达式时

ln( exp(x_1) + exp(x_2) + ... )

你可能会遇到两种问题：

• exp(x_i) 可能会溢出（x_i 太大），导致无法相加的数字
• exp(x_i) 可能会下溢（x_i 太小），导致一堆零

如果所有值都很大或都很小，我们可以除以某个 exp(const) 并在 ln 的外部添加 const 以获得相同的值。因此，如果我们能选择恰当的 const，就可以将值移动到某个范围内以防止溢出/下溢。

提问者的问题是，为什么我们选择 max(x_i) 作为这个常数，而不是任何其他值？为什么我们不递归地执行此计算，从每个子集中选取最大值并重复计算对数？

答案是：因为这不重要。

原因是什么呢？假设 x_1 = 10 是大数，x_2 = -10 是小数。（这些数字甚至在数量级上都不是很大，对吧？）那么表达式

ln( exp(10) + exp(-10) ) 

如果您不相信我，可以试试，ln(exp（x_1）+exp（x_2）+…)通常会给出非常接近于max（x_i） 的值，如果某些特定的x_i比其他所有数都大得多。 （顺便说一句，这种函数形式在渐近意义下实际上让您能够从一组数字中数学选择最大值。）

因此，我们选择最大值而不是任何其他值的原因是较小的值几乎不会影响结果。 如果它们下溢，则它们本来就太小了，无法影响总和，因为它将受到最大数字及其接近值的支配。在计算术语中，小数字的贡献将在计算ln后小于一个ulp。 因此，如果它们最终会丢失，那么没有必要递归地计算较小值的表达式，因为它们将在最终结果中丢失。

如果您想真正挑剔地实现此操作，您可以除以exp（max（x_i）- some_constant）或类似物以“居中”结果值约为1，以避免溢出和下溢，并且可能会在结果中提供几个额外的数字精度。 但是，避免溢出比避免下溢更重要，因为前者确定了结果而后者没有，因此按照这种方式做要简单得多。

- Andrew Mao

Answer 2

递归并不能更好地解决这个问题。问题在于你需要确保有限精度算术不会将噪声淹没掉答案。通过单独处理最大值，可以确保任何垃圾数据在最终答案中都很小，因为其中最重要的组成部分是保证能够得到的。

对于我那啰嗦的解释，我道歉了。您可以自己尝试使用一些数字（一个合理的起始列表可能是[1E-5，1E25，1E-5]），以此来感受它所产生的效果。

- Nicholas Wilson

Answer 3

根据你的定义，你的“recursive”函数将永远不会停止。这是因为“((x-maxlog) for x in listOfLogs)”仍然具有与“listOfLogs”相同数量的元素。

我认为这也不容易修复，除非对性能或精度（与非递归版本相比）产生重大影响。

- RBarryYoung

---

是的，就像我写的那样，它不会结束。而且两次迭代的结果与第一次迭代相同，只是将最大值替换为第二大的值。 - Jacob