Katastrofalne anulowanie w logarytmie

Próbuję zaimplementować następującą funkcję w zmiennoprzecinkowym podwójnej precyzji z niskim błędem względnym :

l o g s u m (x, y) = \log (\exp (x) + \exp (y))

$\mathrm{logsum}(x,y) = \log(\exp(x) + \exp(y))$

Jest to szeroko stosowane w aplikacjach statystycznych w celu dodania prawdopodobieństw lub gęstości prawdopodobieństwa, które są reprezentowane w przestrzeni dziennika. Oczywiście albo $\exp(x)$ albo $\exp(y)$ mogą łatwo przepełnić lub niedopełnić, co byłoby złe, ponieważ przestrzeń dziennika jest używana przede wszystkim do uniknięcia niedopełnienia. To typowe rozwiązanie:

l o g s u m (x, y) = x + l o g 1 p (\exp (y - x))

$\mathrm{logsum}(x,y) = x + \mathrm{log1p}(\exp(y - x))$

Anulowanie z $y-x$ ma miejsce, ale jest złagodzone przez $\exp$ . Zdecydowanie gorsze jest, gdy $x$ i $\mathrm{log1p}(\exp(y - x))$ są blisko. Oto wykres błędu względnego:

wprowadź opis zdjęcia tutaj

Wykres jest odcięty w podkreślać kształt krzywej , o których występuje odwołanie. Widziałem błędu do i podejrzewa, że robi się znacznie gorzej. (FWIW, funkcja „prawdy gruntu” jest implementowana za pomocą pływaków MPFR o dowolnej dokładności z 128-bitową precyzją). $10^{-14}$ $\mathrm{logsum}(x,y) = 0$ $10^{-11}$

Próbowałem innych przeformułowań, wszystkie z tym samym rezultatem. Przy jako zewnętrzna ekspresji sam błąd występuje poprzez rejestr coś w pobliżu 1. Z , co zewnętrzna ekspresji, odwołanie się dzieje w wewnętrznej ekspresji. $\log$ $\mathrm{log1p}$

Teraz absolutnym błędu jest bardzo mały, więc jest bardzo mały błąd względny (szerokość epsilon). Można argumentować, że ponieważ użytkownik jest naprawdę zainteresowany prawdopodobieństwa (prawdopodobieństw nie log), to straszny błąd względny nie jest problemem. Prawdopodobnie zwykle tak nie jest, ale piszę funkcję biblioteczną i chciałbym, aby jej klienci mogli liczyć na błąd względny niewiele gorszy niż błąd zaokrąglania. $\exp(\mathrm{logsum}(x,y))$ $\mathrm{logsum}$

Wygląda na to, że potrzebuję nowego podejścia. Co to może być

floating-point stability numerics

— Neil Toronto
źródło

Nie rozumiem twojego ostatniego akapitu. „w epsilon” nic dla mnie nie znaczy. Masz na myśli jednostkę na ostatnim miejscu ? Jeśli chodzi o użytkowników zainteresowanych prawdopodobieństwem, mały błąd prawdopodobieństwa dziennika spowoduje duży błąd prawdopodobieństwa, więc tak nie jest.

— Aron Ahmadia

Czy z ciekawości próbowałeś wykorzystać to, co najlepsze z dwóch metod i wykreślić błąd? Wtedy wszystko, czego potrzebujesz, to właściwa logika do wykrycia, w którym jesteś przypadku (mam nadzieję, że i tak będzie mniej kosztowna lub część wymaganego kosztu algorytmu), a następnie przejście na odpowiednią metodę.

— Aron Ahmadia

@AronAhmadia: „W epsilonie” oznacza błąd względny mniejszy niż epsilon zmiennoprzecinkowy podwójnej precyzji, który wynosi około 2,22e-16. Dla normalnych (tj. Nie nienormalnych) pływaków odpowiada to około ulp. Ponadto, jeśli

jest błędem bezwzględnym

, to błąd względny

wynosi

, co jest prawie funkcją tożsamości zbliżoną do zera. IOW, mały błąd bezwzględny dla

implikuje mały błąd względny dla

a

$a$

x

$x$

\exp (x)

$\exp(x)$

\exp (a) - 1

$\exp(a)-1$

x

$x$

\exp (x)

$\exp(x)$

— Neil Toronto,

Dodatek: Gdy błąd bezwzględny

jest bliski zeru. Kiedy

, na przykład, masz rację: względne eksploduje.

a

$a$

a > 1

$a > 1$

— Neil Toronto,

Formuła powinna być stabilna numerycznie. Uogólnia się na obliczenie stabilne numerycznie

l o sol s u m (x, y) = max (x, y) + l o sol 1 p (\exp (- abs (x - y))

$\mathrm{logsum}(x,y)=\max(x,y)+\mathrm{log1p}(\exp(-\operatorname{abs}(x-y))$

\log \sum_{ja} {mi}^{x_{ja}} = ξ + \log \sum_{ja} {mi}^{x_{ja} - ξ}, ξ = max_{ja} x_{ja}

$\log \sum_i e^{x_i} = \xi+ \log\sum_i e^{x_i-\xi},~~~\xi=\max_i x_i$

W przypadku, gdy suma logarytmiczna jest bardzo bliska zeru i chcesz wysokiej dokładności względnej, prawdopodobnie możesz użyć przy użyciu dokładna (tj. ponad dwukrotna precyzja) implementacja

l o sol s u m (x, y) = max (x, y) + l mi x p (x - y)

$\mathrm{logsum}(x,y)=\max(x,y)+\mathrm{lexp}(x-y)$

l mi x p (z) : = \log (1 + {mi}^{- | z |})

$\mathrm{lexp}(z):=\log(1+e^{-|z|})$ co jest prawie liniowe dla małego

z

$z$

— Arnold Neumaier
źródło

Pod względem błędu bezwzględnego tak jest. Pod względem błędu względnego jest okropnie, gdy wyjście jest bliskie zeru.

— Neil Toronto,

x

$x$

y

$y$

Dla x = -0,775 i y = -0,6175, otrzymuję błąd 62271 ulps i błąd względny 1,007e-11.

— Neil Toronto,

Oblicz bardzo dokładne punkty danych w zakresie zainteresowania - co najmniej dwa różne zakresy są potrzebne z powodu zachowania asymptotycznego. Można użyć wyrażenia definiującego dla z nie zbliżonego do zera. Dla wyjątkowego zakresu dopasuj racjonalną funkcję o wystarczająco wysokim stopniu, aby uzyskać pożądaną dokładność. Aby uzyskać stabilność numeryczną, użyj wielomianów Bernsteina lub wielomianów Tchebycheva w liczniku i mianowniku, dostosowanych do przedziału zainteresowania. Na koniec rozszerz się na ciągłą frakcję i dowiedz się, jak bardzo można skrócić współczynniki bez utraty dokładności.

— Arnold Neumaier,

l = l (z)

$l=l(z)$

m

$m$