Как реализовать численно стабильное взвешенное логарифмическое выражение?

Какой наиболее численно устойчивый способ расчета:

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

Как мне обобщить это на взвешенную версию?


math language-agnostic numerical-stability
person Neil G    schedule 12.07.2015    source источник

Ответы (2)


arrow_upward
5
arrow_downward

Вы можете использовать исходную функцию logaddexp для этой цели, если вы перепишете взвешенное выражение как,

новое выражение logadd

Это эквивалентно,

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

который должен быть таким же численно стабильным, как и исходная реализация logaddexp.

Примечание. я имею в виду numpy.logaddexp, которая принимает x и y, а не x и exp_y, как вы упомянули в вопросе.

person rth    schedule 12.07.2015
comment
Похоже, это, вероятно, лучше, чем то, что я сделал, что я добавлю в качестве ответа для сравнения. - person Neil G; 12.07.2015
comment
Для тех, кто читает, я проверил это с помощью библиотеки произвольной точности mpmath и обнаружил, что это намного лучше, чем мое решение. - person Neil G; 16.07.2015
comment
@NeilG Да, я полагаю, что независимо от того, как вы его перепишете, вы все равно потеряете точность / переполнение с 64-битными числами с плавающей запятой и т. д., когда берете экспоненту больших чисел и вычисляете журнал. mpmath кажется здесь хорошим выбором, хотя и будет медленнее. - person rth; 16.07.2015
comment
Я не использую mpmath в своем коде. Я просто использовал его для измерения средней ошибки каждого из наших решений, чтобы знать, какое из них на самом деле лучше. mpmath дает ответ произвольной точности. Моя интуиция при написании моего решения заключалась в том, чтобы использовать log1p как можно больше. Это оказалось намного хуже. - person Neil G; 16.07.2015

arrow_upward
1
arrow_downward 

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))
person Neil G    schedule 12.07.2015