如何在不溢出的情况下计算log(x - y)，已知log(x)和log(y)？

Question

以下函数计算log(x + y)，给定log(x)和log(y)的值，如果x或y非常大或非常小，则避免溢出或下溢：

double log_add(double logx, double logy)
{
    return max(logx, logy) + log1p(exp(-fabs(logx - logy)));
}

必须有一个类似的log_sub函数来计算log(x - y)。它是什么？

更一般地，我需要计算log(x - y - z)，已知log(x)、log(y)和log(z)。通过log_add和log_sub，我可以在两个步骤中计算log(x - y - z)，但也许有一种更优化的方法吗？

- a06e

---

@LưuVĩnhPhúc 不对。你否定的是y的对数，而不是y本身。 - a06e

抱歉，我误读了，这些值是日志。 - phuclv

Answer 1

为什么不直接使用对数恒等式呢：

double log_add(double logx, double logy) {
    return logx + log1p(exp(logy - logx));
}

double log_sub(double logx, double logy) {
    return logx + log1p(-exp(logy - logx));
}

针对您的具体情况：

// log(x - y - z) given the three logs
double log_xyz(double logx, double logy, double logz) {
    return logx + log1p(-exp(logy - logx) - exp(logz - logx));
}

- Barry

---

由于存在对 exp() 的调用，我认为该代码运行时存在实际的溢出风险。例如，当 logy = 708，logx = -3 时，在 IEEE-754 双精度算术中会导致溢出。 - njuffa

@njuffa 我只是提供了身份信息 :) 这就是它们。 - Barry

请参阅 http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf 了解有关log1pexp的详细信息。 - Simon Byrne

很不幸，CRAN目前似乎已经宕机，请尝试点击此链接http://cran.stat.ucla.edu/web/packages/Rmpfr/vignettes/log1mexp-note.pdf。 - Simon Byrne

@njuffa 我假设这段代码的意图是在logx > logy的情况下被调用，但是在这种情况下我们有很大的机会遇到下溢的问题 :) - aka.nice

如果您在所有函数中说明x、y、z的相对大小或logx、logy、logz的大小关系，那将会很有帮助。我已经有点相信，存在一些假设使得您的答案是正确和有用的，但很难确定。 - Don Hatch