如何实现一个数值稳定的加权对数求和？（涉及IT技术）

Question

什么是最稳定的计算方式来计算：

log[(wx * exp(x) + wy * exp_y)/(wx + wy)]

这里的权重 wx、wy > 0 是什么意思？

如果没有权重，这个函数就是 logaddexp，可以使用NumPy在Python中实现：

tmp = x - y
return np.where(tmp > 0,
                x + np.log1p(np.exp(-tmp)),
                y + np.log1p(np.exp(tmp)))

我该如何将它推广到加权版本？

- Neil G

---

scipy.special.logsumexp具有所需的行为，其中包括一个权重参数b。 - Wouter

Answer 1

您可以使用原始的logaddexp函数来实现此目的，如果您将加权表达式重写为：

这等价于：

logaddexp( x + log(w_x), y + log(w_y) ) - log(w_x + w_y)

应该保持与原始的logaddexp实现一样的数值稳定性。

注意：我指的是接受x和y作为输入参数的numpy.logaddexp函数，而不是你在问题中提到的x和exp_y。

- rth

---

看起来这可能比我做的更好，我会将其作为答案添加以供比较。 - Neil G

对于任何读者，我使用了任意精度库“mpmath”进行测试，并发现它比我的解决方案要好得多。 - Neil G

@NeilG 是的，我想无论你如何重写它，在计算大数的指数和对数时，使用64位浮点数等仍然会失去精度/溢出。mpmath 在这里似乎是一个不错的选择，尽管速度会慢一些。 - rth

我在我的代码中没有使用mpmath。我只是用它来测量每个解决方案的平均误差，以便我知道哪个实际上更好。mpmath提供了任意精度的答案。当我编写我的解决方案时，我的直觉是尽可能多地使用log1p。结果证明这样做要糟糕得多。 - Neil G

Answer 2

def weighted_logaddexp(x, wx, y, wy):
    # Returns:
    #   log[(wx * exp(x) + wy * exp_y)/(wx + wy)]
    #   = log(wx/(wx+wy)) + x + log(1 + exp(y - x + log(wy)-log(wx)))
    #   = log1p(-wy/(wx+wy)) + x + log1p((wy exp_y) / (wx exp(x)))
    if wx == 0.0:
        return y
    if wy == 0.0:
        return x
    total_w = wx + wy
    first_term = np.where(wx > wy,
                          np.log1p(-wy / total_w),
                          np.log1p(-wx / total_w))
    exp_x = np.exp(x)
    exp_y = np.exp(y)
    wx_exp_x = wx * exp_x
    wy_exp_y = wy * exp_y
    return np.where(wy_exp_y < wx_exp_x,
                    x + np.log1p(wy_exp_y / wx_exp_x),
                    y + np.log1p(wx_exp_x / wy_exp_y)) + first_term

---

这是我比较这两个解决方案的方法：

import math
import numpy as np
import mpmath as mp
from tools.numpy import weighted_logaddexp

def average_error(ideal_function, test_function, n_args):
    x_y = [np.linspace(0.1, 3, 20) for _ in range(n_args)]
    xs_ys = np.meshgrid(*x_y)

    def e(*args):
        return ideal_function(*args) - test_function(*args)
    e = np.frompyfunc(e, n_args, 1)
    error = e(*xs_ys) ** 2
    return np.mean(error)


def ideal_function(x, wx, y, wy):
    return mp.log((mp.exp(x) * wx + mp.exp(y) * wy) / mp.fadd(wx, wy))

def test_function(x, wx, y, wy):
    return np.logaddexp(x + math.log(wx), y + math.log(wy)) - math.log(wx + wy)

mp.prec = 100
print(average_error(ideal_function, weighted_logaddexp, 4))
print(average_error(ideal_function, test_function, 4))

- Neil G