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

Question

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

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)

Но я все еще получаю "inf", чего я хочу избежать.

Я пытался обрезать его, используя 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) 

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

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

Answer 1

Только значительно расширенные пакеты или пакеты произвольной точности смогут справиться с огромными различиями в количестве. Экспонента самых больших и самых отрицательных чисел в 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

Answer 2

Если проблемы потери цифр в окончательных результатах из-за сильно различающихся порядков величин добавленных терминов не являются проблемой, можно также математически преобразовать 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))) численно равен нулю, так как экспонента наибольшего числа здесь полностью доминирует в сумме.

Answer 3

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