Comment ajouter de gros termes exponentiels de manière fiable sans erreurs de débordement?

Un problème très courant dans Markov Chain Monte Carlo implique le calcul de probabilités qui sont la somme de grands termes exponentiels,

$e^{a_1} + e^{a_2} + ...$

où les composants d' peut varier de très petit à très grand. Mon approche a été de factoriser le plus grand terme exponentiel pour que:K : = max i ( a i$a$$K := \max_{i}(a_{i})$

e un ' ≡ e a 1 + e a 2 + . .

Cette approche est raisonnable si tous les éléments d' $a$ sont grands, mais ce n'est pas une bonne idée s'ils ne le sont pas. Bien sûr, les plus petits éléments ne contribuent pas de toute façon à la somme en virgule flottante, mais je ne sais pas comment les gérer de manière fiable. En code R, mon approche ressemble à:

if ( max(abs(a)) > max(a) )
  K <-  min(a)
else
  K <- max(a)
ans <- log(sum(exp(a-K))) + K

Cela semble un problème assez courant qu'il devrait y avoir une solution standard, mais je ne suis pas sûr de ce que c'est. Merci pour toutes suggestions.

Il existe une solution simple avec seulement deux passages dans les données:

Calculez d'abord

ce qui vous dit que, s'il y a termes, alors $n$

Comme vous n'avez probablement pas aussi grand que , vous ne devriez pas vous soucier de déborder dans le calcul de en double précision .$n$$10^{20}$

Ainsi, calculez et votre solution est .$\tau$$e^K \tau$

Pour conserver la précision pendant que vous ajoutez des doubles ensemble, vous devez utiliser Kahan Summation , c'est l'équivalent logiciel d'avoir un registre de retenue.

C'est très bien pour la plupart des valeurs, mais si vous obtenez un débordement, vous atteignez les limites de la double précision IEEE 754 qui serait d'environ . À ce stade, vous avez besoin d'une nouvelle représentation. Vous pouvez détecter un débordement au moment de l'ajout et également détecter des exposants trop grands pour être évalués par . À ce stade, vous pouvez modifier l'interprétation d'un double en décalant l'exposant et en gardant une trace de ce décalage.$e^{709.783}$doubleMax - sumSoFar < valueToAddexponent > 709.783

Ceci est pour la plupart similaire à votre approche de décalage d'exposant, mais cette version est conservée en base 2 et ne nécessite pas de recherche initiale pour trouver le plus grand exposant. D'où .$\mathit{value}\times2^\mathit{shift}$

#!/usr/bin/env python
from math import exp, log, ceil

doubleMAX = (1.0 + (1.0 - (2 ** -52))) * (2 ** (2 ** 10 - 1))

def KahanSumExp(expvalues):
  expvalues.sort() # gives precision improvement in certain cases 
  shift = 0 
  esum = 0.0 
  carry = 0.0 
  for exponent in expvalues:
    if exponent - shift * log(2) > 709.783:
      n = ceil((exponent - shift * log(2) - 709.783)/log(2))
      shift += n
      carry /= 2*n
      esum /= 2*n
    elif exponent - shift * log(2) < -708.396:
      n = floor((exponent - shift * log(2) - -708.396)/log(2))
      shift += n
      carry *= 2*n
      esum *= 2*n
    exponent -= shift * log(2)
    value = exp(exponent) - carry 
    if doubleMAX - esum < value:
      shift += 1
      esum /= 2
      value /= 2
    tmp = esum + value 
    carry = (tmp - esum) - value 
    esum = tmp
  return esum, shift

values = [10, 37, 34, 0.1, 0.0004, 34, 37.1, 37.2, 36.9, 709, 710, 711]
value, shift = KahanSumExp(values)
print "{0} x 2^{1}".format(value, shift)

Votre approche est solide.

$K$$K$

Il existe un package R qui fournit une implémentation rapide et efficace de "l'astuce log-sum-exp"

http://www.inside-r.org/packages/cran/matrixStats/docs/logSumExp

La fonction logSumExp accepte un vecteur numérique lX et génère un journal (sum (exp (lX))) tout en évitant les problèmes de débordement et de débordement en utilisant la méthode que vous avez décrite.