Переполнение в numpy.exp()

Мне нужно вычислить экспоненту следующего массива для моего проекта:

w  = [-1.52820754859, -0.000234000845064, -0.00527938881237, 5797.19232191, -6.64682108484,
       18924.7087966, -69.308158911, 1.1158892974, 1.04454511882, 116.795573742]

Но у меня было переполнение из-за номера 18924.7087966.

Цель состоит в том, чтобы избежать использования дополнительных пакетов, таких как bigfloat (кроме «numpy»), и получить близкий результат (который имеет небольшую относительную ошибку).

1. До сих пор я пытался использовать более высокую точность (например, float128):

def getlogZ_robust(w):

    Z = sum(np.exp(np.dot(x,w).astype(np.float128)) for x in iter_all_observations())
    return np.log(Z)

Но я все еще получаю «информацию», чего я хочу избежать.

  1. Я попытался обрезать его с помощью nump.clip():

    def getlogZ_robust(w):

    Z = sum(np.exp(np.clip(np.dot(x,w).astype(np.float128),-11000, 11000)) for x in iter_all_observations())
    return np.log(Z)

Но относительная ошибка слишком велика.

Можете ли вы помочь мне решить эту проблему, если это возможно?


python numpy overflow exp exponential
person erenerogullari    schedule 10.11.2019    source источник
comment
Можете ли вы использовать SciPy? Если это так, проверьте logsumexp. Если нет, вы можете получить исходный код с github.com/ scipy/scipy/blob/master/scipy/special/_logsumexp.py или немного поискать здесь, в stackoverflow. Вероятно, существует множество реализаций одного и того же метода, разбросанных по ответам stackoverflow.   -  person Warren Weckesser    schedule 10.11.2019
comment
Я не могу использовать SciPy, но я проверил logsumexp и обнаружил, что у numpy также есть аналогичная функция под названием logaddexp, которая вычисляет logaddexp(x1, x2) == log(exp(x1) + exp(x2)) без явного вычисления промежуточных значений exp(). Таким образом, он избегает переполнения. Так что все-таки спасибо.   -  person erenerogullari    schedule 10.11.2019

Ответы (3)


arrow_upward
4
arrow_downward

Только значительно расширенные пакеты или пакеты произвольной точности смогут справиться с огромными различиями в числах. Экспонента самого большого и самого отрицательного числа в w отличается на 8000 (!) порядков. float (то есть двойная точность) имеет «всего» 15 цифр точности (что означает, что 1+1e-16 численно равно 1), так что добавление маленьких чисел к огромной экспоненте наибольшего числа не имеет никакого эффекта. На самом деле, exp(18924.7087966) настолько велико, что доминирует в сумме. Ниже приведен скрипт, выполняющий сумму с повышенной точностью в mpmath: отношение суммы экспонент и exp(18924.7087966) в основном равно 1.

w  = [-1.52820754859, -0.000234000845064, -0.00527938881237, 5797.19232191, -6.64682108484,
       18924.7087966, -69.308158911, 1.1158892974, 1.04454511882, 116.795573742]

u = min(w)
v = max(w)

import mpmath
#using plenty of precision
mpmath.mp.dps = 32768
print('%.5e' % mpmath.log10(mpmath.exp(v)/mpmath.exp(u)))
#exp(w) differs by 8000 orders of magnitude for largest and smallest number

s = sum([mpmath.exp(mpmath.mpf(x)) for x in w])

print('%.5e' % (mpmath.exp(v)/s))
#largest exp(w) dominates such that ratio over the sums of exp(w) and exp(max(w)) is approx. 1
person Gianluca    schedule 10.11.2019

arrow_upward
3
arrow_downward

Если проблемы с потерей цифр в окончательных результатах из-за очень разных порядков величин добавленных членов не вызывают беспокойства, можно также математически преобразовать log сумм по экспонентам следующим образом, избегая exp больших чисел:

log(sum(exp(w)))
= log(sum(exp(w-wmax)*exp(wmax)))
= wmax + log(sum(exp(w-wmax)))

В питоне:

import numpy as np
v = np.array(w)
m = np.max(v)
print(m + np.log(np.sum(np.exp(v-m))))

Обратите внимание, что np.log(np.sum(np.exp(v-m))) численно равен нулю, так как экспонента наибольшего числа полностью доминирует над суммой здесь.

person Gianluca    schedule 10.11.2019

arrow_upward
1
arrow_downward

В Numpy есть функция logaddexp, которая вычисляет

logaddexp(x1, x2) == log(exp(x1) + exp(x2))

без явного вычисления промежуточных значений exp(). Таким образом, он избегает переполнения. Итак, вот решение:

def getlogZ_robust(w):

    Z = 0
    for x in iter_all_observations():
        Z = np.logaddexp(Z, np.dot(x,w))
    return Z
person erenerogullari    schedule 10.11.2019