How to implement numerically stable weighted logaddexp?

Question

What is the most numerically stable way to calculate:

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

where are the weights wx, wy > 0 ?

Without weights, this function is logaddexp and can be implemented in Python with NumPy as:

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

How do I generalize this to a weighted version?

Answer 1

You can use the original logaddexp function for this purpose if you rewrite the weighted expression as,

It is equivalent

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

which should be as numerically robust as the original implementation of logaddexp .

Note: I am referring to the function numpy.logaddexp , which accepts x and y , not x and exp_y , as you mentioned in the question.

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 

---

Here is how I compared the two solutions:

 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))