Избегайте недополнения, используя exp и минимальное положительное значение float128 в numpy

Я пытаюсь рассчитать следующее соотношение: w(i) / (sum(w(j)), где w обновляются с использованием функции экспоненциального убывания, т.е. w(i) = w(i) * exp(-k), k являются положительными параметрами. Все числа неотрицательны. Затем это отношение используется в формуле (умножить на константу и добавить другую константу). Как и ожидалось, вскоре я столкнулся с проблемами недополнения.

Я думаю, это происходит часто, но может ли кто-нибудь дать мне несколько ссылок о том, как с этим бороться? Я не нашел подходящего преобразования, поэтому одну вещь, которую я попытался сделать, это установить минимальное положительное число в качестве порога безопасности, но мне не удалось найти минимальное положительное число с плавающей запятой (я представляю числа в numpy.float128). Как я могу получить минимальное положительное такое число на моей машине? Код выглядит следующим образом:

w = np.ones(n, dtype='float128')
lt = np.ones(n)
for t in range(T):
    p = (1-k) * w / w.sum() + (k/n)
    # Process a subset of the n elements, call it set I, j is some range()
    for i in I: 
        s = p[list(j[i])].sum()
        lt /= s
        w[s] *= np.exp(-k * lt)

где k — некоторая константа в (0,1), а n — длина массива


python numpy numerical-methods underflow numerical-analysis
person user90772    schedule 30.10.2015    source источник
comment
numpy.float128 — очень вводящее в заблуждение имя — его реальная точность такая же, как и у вашего компилятора C, который называется long double, что составляет 80 бит на x86 (для ясности вместо этого рекомендуется использовать псевдоним numpy.longdouble). Чтобы получить наименьшее представимое положительное число с плавающей запятой, вы можете использовать np.finfo(numpy.float128).tiny который на моей машине равен 3.3621031431120935063e-4932, но более полезным порогом, вероятно, будет машинный эпсилон (eps).   -  person ali_m    schedule 30.10.2015
comment
Не могли бы вы показать нам свой код? Я подозреваю, что, вероятно, есть лучший способ вычисления этого отношения.   -  person ali_m    schedule 30.10.2015
comment
Спасибо @ali_m. Я добавил фрагмент кода, я не мог добавить его полностью, но идея в том, что я обрабатываю другое подмножество ws во времени и уменьшаю их, используя np.exp().   -  person user90772    schedule 02.11.2015

Ответы (1)


arrow_upward
3
arrow_downward

При работе с экспоненциально малыми числами обычно лучше работать в пространстве журнала. Например, log(w*exp(-k)) = log(w) - k, у которого не будет проблем с переполнением/недостаточным значением, если только k не является экспоненциально большим или w не равно нулю. И если w равно нулю, numpy правильно вернет -inf. Затем, выполняя сумму, вы выносите за скобки наибольший член:

log_w = np.log(w) - k
max_log_w = np.max(log_w)
# Individual terms in the following may underflow, but then they wouldn't
# contribute to the sum anyways.
log_sum_w = max_log_w + np.log(np.sum(np.exp(log_w - max_log_w)))
log_ratio = log_w - log_sum_w

Это, вероятно, не совсем то, что вам нужно, поскольку вы можете просто полностью исключить k (при условии, что это константа, а не массив), но это должно помочь вам.

Scikit-learn реализует нечто подобное с extmath.logsumexp, но в основном это то же самое, что и выше.

person clwainwright    schedule 30.10.2015
comment
Также есть numpy.logaddexp и < a href="http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.misc.logsumexp.html" rel="nofollow noreferrer">scipy.misc.logsumexp - person ali_m; 31.10.2015
comment
Не знал о scipy.misc.logsumexp, похоже, это было бы особенно полезно для ОП. Спасибо! - person clwainwright; 31.10.2015
comment
np.logaddexp также является ufunc, поэтому вы можете использовать np.logaddexp.reduce для суммирования по оси массива. - person ali_m; 31.10.2015