如何实现数值稳定的加权 logaddexp?
How do I implement a numerically stable weighted logaddexp?
什么是数值最稳定的计算方式:
log[(wx * exp(x) + wy * exp_y)/(wx + wy)]
权重在哪里 wx, wy > 0
?
没有权重,此函数是 logaddexp
并且可以在 Python 中使用 NumPy 实现为:
tmp = x - y
return np.where(tmp > 0,
x + np.log1p(np.exp(-tmp)),
y + np.log1p(np.exp(tmp)))
我应该如何将其推广到加权版本?
如果将加权表达式重写为
,则可以将原始 logaddexp
函数用于此目的
这相当于,
logaddexp( x + log(w_x), y + log(w_y) ) - log(w_x + w_y)
这应该与原始 logaddexp
实施一样在数值上稳定。
注意: 我指的是接受 x
和 y
的 numpy.logaddexp
函数,而不是 x
和 exp_y
,正如您在问题中提到的。
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))
什么是数值最稳定的计算方式:
log[(wx * exp(x) + wy * exp_y)/(wx + wy)]
权重在哪里 wx, wy > 0
?
没有权重,此函数是 logaddexp
并且可以在 Python 中使用 NumPy 实现为:
tmp = x - y
return np.where(tmp > 0,
x + np.log1p(np.exp(-tmp)),
y + np.log1p(np.exp(tmp)))
我应该如何将其推广到加权版本?
如果将加权表达式重写为
,则可以将原始logaddexp
函数用于此目的
这相当于,
logaddexp( x + log(w_x), y + log(w_y) ) - log(w_x + w_y)
这应该与原始 logaddexp
实施一样在数值上稳定。
注意: 我指的是接受 x
和 y
的 numpy.logaddexp
函数,而不是 x
和 exp_y
,正如您在问题中提到的。
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))