Wie implementiere ich eine numerisch stabile gewichtete logaddexp?

Question

Was ist die numerisch stabile Art und Weise der Berechnung:

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

wo die Gewichte wx, wy > 0?

Ohne die Gewichte, diese Funktion logaddexp und könnte wie in Python mit NumPy implementiert werden:

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

Wie soll ich verallgemeinern dies auf die gewichtete Version?

Answer 1

Sie konnten die ursprünglichen logaddexp Funktion für so Zweck verwenden, wenn Sie den gewichteten Ausdruck als umschreiben,

Dies äquivalent ist,

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

, die als numerisch stabil sein sollte als die ursprüngliche logaddexp Implementierung.

Hinweis: Ich beziehe mich auf die numpy.logaddexp-Funktion, die in x und y nimmt, nicht x und exp_y, wie Sie in der Frage erwähnen.

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 

---

Hier ist, wie ich die beiden Lösungen miteinander verglichen:

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