Jak niezawodnie dodawać duże wyrażenia wykładnicze bez błędów przepełnienia?

24

Bardzo częstym problemem w Markov Chain Monte Carlo jest obliczanie prawdopodobieństw, które są sumą dużych terminów wykładniczych,

ea1+ea2+...

gdzie składniki może wahać się od bardzo małych do bardzo dużych. Moje podejście na uwzględnieniu największego wykładniczego terminu , aby:aK:=maxi(ai)

e 'e się 1 + wiadomości e się 2 + . . .

a=K+log(ea1K+ea2K+...)
eaea1+ea2+...

Takie podejście jest uzasadnione, jeżeli wszystkie elementy są duże, ale nie taki dobry pomysł, jeśli nie są. Oczywiście, mniejsze elementy i tak nie przyczyniają się do sumy zmiennoprzecinkowej, ale nie jestem pewien, jak sobie z nimi poradzić. W kodzie R moje podejście wygląda następująco:a

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

Wydaje się dość powszechnym problemem, że powinno istnieć standardowe rozwiązanie, ale nie jestem pewien, co to jest. Dziękuję za wszelkie sugestie.

cboettig
źródło
1
To jest rzecz. Google dla „logsumexp”.

Odpowiedzi:

15

Istnieje proste rozwiązanie z dwoma przejściami danych:

Najpierw oblicz

K:=maxiai,

co mówi ci, że jeśli nie ma terminów, to Σ i e in e K .n

ieaineK.

Ponieważ prawdopodobnie nie masz tak blisko jak nawet , nie powinieneś się martwić przepełnieniem obliczeń z podwójną precyzją .10 20 τ : = i e a i - Knn1020

τ:=ieaiKn

Zatem oblicz / a wtedy twoim rozwiązaniem jest .e K ττeKτ

Jack Poulson
źródło
Dzięki za wyraźną notację - ale uważam, że to jest właśnie to, co zaproponowałem (?) Jeśli muszę unikać błędów niedomiaru, gdy niektóre są małe, rozumiem, że potrzebuję podejścia sumowania Kahan zaproponowanego przez @gareth ? ai
cboettig
Ach, teraz widzę, o co ci chodziło. W rzeczywistości nie musisz martwić się o niedomiar, ponieważ dodanie wyjątkowo małych wyników do rozwiązania nie powinno tego zmienić. Jeśli było ich wyjątkowo dużo, najpierw należy zsumować małe wartości.
Jack Poulson
Do downvoter: czy mógłbyś mi powiedzieć, co jest nie tak z moją odpowiedzią?
Jack Poulson
co jeśli masz wiele bardzo małych terminów? Może się zdarzyć, że dla nich. Jeśli istnieje wiele takich terminów, wystąpiłby duży błąd. eaiK0
becko
10

Aby zachować precyzję podczas sumowania podwójnego, musisz użyć Kahan Summation , jest to oprogramowanie równoważne z rejestrem przenoszenia.

Jest to w porządku dla większości wartości, ale jeśli dostajesz przepełnienie, to osiągasz limit podwójnej precyzji IEEE 754, co byłoby około . W tym momencie potrzebujesz nowej reprezentacji. Możesz wykryć przepełnienie w czasie dodawania, a także wykryć duże wykładniki do oceny . W tym momencie możesz zmodyfikować interpretację podwójnego, przesuwając wykładnik i śledząc tę ​​zmianę.e709.783doubleMax - sumSoFar < valueToAddexponent > 709.783

To w przeważającej części jest podobne do twojego podejścia wykładniczego, ale ta wersja jest przechowywana w bazie 2 i nie wymaga wstępnego wyszukiwania w celu znalezienia największego wykładnika. Stąd .value×2shift

#!/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)
Gareth A. Lloyd
źródło
Sumowanie Kahana jest tylko jedną z rodziny metod „kompensowanego sumowania”. Jeśli z jakiegoś powodu Kahan nie działa całkiem dobrze, istnieje wiele innych metod poprawnego poprawiania warunków o różnej wielkości i przeciwnych znakach.
JM
@JM, czy możesz podać mi nazwy tych innych metod, chętnie bym je przeczytał. Dzięki.
Gareth A. Lloyd
1

Twoje podejście jest solidne.

Nie musisz dokładnie znać , wystarczy, aby uniknąć przepełnienia. Być może będziesz w stanie oszacować analitycznie przed wykonaniem próbkowania MCMC.KKK

John D. Cook
źródło
0

Istnieje pakiet R, który zapewnia szybką i wydajną implementację „sztuczki log-sum-exp”

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

Funkcja logSumExp akceptuje wektor numeryczny lX i log danych wyjściowych (suma (exp (lX))), unikając problemów z przepełnieniem i przepełnieniem przy użyciu opisanej metody.

kantai
źródło