Katastrofalne anulowanie w logarytmie

18

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

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:

logsum(x,y)=x+log1p(exp(yx))

Anulowanie z yx ma miejsce, ale jest złagodzone przez exp . Zdecydowanie gorsze jest, gdy x i log1p(exp(yx)) są blisko. Oto wykres błędu względnego:

wprowadź opis zdjęcia tutaj

Wykres jest odcięty w podkreślać kształt krzywej l o, g y U m ( x , y ) = 0 , o których występuje odwołanie. Widziałem błędu do 10 - 11 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ą).1014logsum(x,y)=01011

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 l O g 1 s , co zewnętrzna ekspresji, odwołanie się dzieje w wewnętrznej ekspresji.loglog1p

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 l o g s u m 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(logsum(x,y))logsum

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

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 x , to błąd względny exp ( x ) wynosi exp ( a ) - 1 , co jest prawie funkcją tożsamości zbliżoną do zera. IOW, mały błąd bezwzględny dla x implikuje mały błąd względny dla exp ( x ) . zaxexp(x)exp(za)-1xexp(x)
Neil Toronto,
Dodatek: Gdy błąd bezwzględny jest bliski zeru. Kiedy > 1 , na przykład, masz rację: względne eksploduje. zaza>1
Neil Toronto,

Odpowiedzi:

12

Formuła powinna być stabilna numerycznie. Uogólnia się na obliczenie stabilne numerycznie log i e x i = ξ + log i e x i

losolsum(x,y)=max(x,y)+losol1p(exp(-abs(x-y))
logjamixja=ξ+logjamixja-ξ,   ξ=maxjaxja

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 e x p ( z ) : = log ( 1 + e - | z |

losolsum(x,y)=max(x,y)+lmixp(x-y)
lmixp(z): =log(1+mi-|z|)
co jest prawie liniowe dla małego .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,
xy
Dla x = -0,775 i y = -0,6175, otrzymuję błąd 62271 ulps i błąd względny 1,007e-11.
Neil Toronto,
1
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,
1
l=l(z)m