Kiedy należy używać log1p i expm1?

30

Mam proste pytanie, które jest naprawdę trudne dla Google (oprócz kanonicznego Co każdy informatyk powinien wiedzieć o arytmetyki zmiennoprzecinkowej ).

Kiedy należy używać funkcji takich jak log1plub expm1należy używać zamiast logi exp? Kiedy nie należy ich używać? Czym różnią się różne implementacje tych funkcji pod względem ich wykorzystania?

Tim
źródło
2
Witamy w Scicomp.SE! To bardzo rozsądne pytanie, ale łatwiej byłoby na nie odpowiedzieć, gdybyś wyjaśnił trochę, o czym log1p mówisz (szczególnie, jak to jest zaimplementowane, więc nie musimy zgadywać).
Christian Clason,
4
W przypadku argumentów o wartościach rzeczywistych należy stosować log1p i expm1 ( x ), gdy x jest małe, np. Gdy 1 + x = 1 z dokładnością zmiennoprzecinkową. Zobacz np. Docs.scipy.org/doc/numpy/reference/generated/numpy.expm1.html i docs.scipy.org/doc/numpy/reference/generated/numpy.log1p.html . (x)(x)x1+x=1
GoHokies,
@ChristianClason dzięki, mam na myśli głównie C ++ std lub R, ale kiedy pytasz, zaczynam myśleć, że poznanie różnic we wdrożeniach byłoby również bardzo interesujące.
Tim
Oto przykład: scicomp.stackexchange.com/questions/8371/...
Juan M. Bello-Rivas
1
@ użytkownik2186862 „gdy jest małe” jest poprawne, ale nie tylko „gdy 1 + x = 1 z dokładnością zmiennoprzecinkową” (co zdarza się dla x 10 - 16 w zwykłej arytmetyce podwójnej precyzji). Stron dokumentacji, którą powiązana pokazać, że są one przydatne już dla x 10 - 10 , na przykład. x1+x=1x1016x1010
Federico Poloni,

Odpowiedzi:

25

Wszyscy wiemy, że

exp(x)=n=0xnn!=1+x+12x2+
oznacza, że ​​dla|x|1, mamyexp(x)1+x. Oznacza to, że jeśli musimy oceniać w zmiennoprzecinkowymexp(x)1, dla|x|1może wystąpić katastrofalne anulowanie.

Można to łatwo zademonstrować w Pythonie:

>>> from math import (exp, expm1)

>>> x = 1e-8
>>> exp(x) - 1
9.99999993922529e-09
>>> expm1(x)
1.0000000050000001e-08

>>> x = 1e-22
>>> exp(x) - 1
0.0
>>> expm1(x)
1e-22

Dokładne wartości

exp(108)1=0.000000010000000050000000166666667083333334166666668exp(1022)1=0.000000000000000000000100000000000000000000005000000

Zasadniczo „dokładna” implementacja expi expm1powinna być poprawna do nie więcej niż 1ULP (tj. Jedna jednostka ostatniego miejsca). Ponieważ jednak osiągnięcie tej dokładności powoduje powstanie „wolnego” kodu, czasami dostępna jest szybka, mniej dokładna implementacja. Na przykład w CUDA mamy expfi expm1f, gdzie foznacza szybko. Według przewodnika programowania CUDA C, aplikacja. Dexpf ma błędu 2ULP.

Jeśli nie przejmujesz się błędami w kolejności kilku ULPS, zwykle różne implementacje funkcji wykładniczej są równoważne, ale uważaj, że błędy mogą być gdzieś ukryte ... (Pamiętasz błąd Pentium FDIV ?)

Jest więc całkiem jasne, że expm1należy użyć do obliczenia exp(x)1 dla małego x . Używanie go do ogólnego x nie jest szkodliwe, ponieważ expm1oczekuje się, że będzie dokładny w całym zakresie:

>>> exp(200)-1 == exp(200) == expm1(200)
True

(W powyższym przykładzie 1 jest znacznie poniżej 1ULP exp(200) , więc wszystkie trzy wyrażenia zwracają dokładnie tę samą liczbę zmiennoprzecinkową.)

Podobna dyskusja dotyczy funkcji odwrotnych logi log1pponieważ log(1+x)x dla |x|1 .

Stefano M.
źródło
1
Ta odpowiedź była już zawarta w komentarzach do pytania PO. Czułem jednak, że warto podać dłuższe (choć podstawowe) konto tylko dla jasności, w nadziei, że przyda się niektórym czytelnikom.
Stefano M
OK, ale wtedy można po prostu stwierdzić „więc zawsze mogę użyć expm1 zamiast exp” ...
Tim
1
@ time, twój wniosek jest błędny: zawsze możesz użyć expm1(x)zamiast exp(x)-1. Oczywiście exp(x) == exp(x) - 1nie ma to miejsca.
Stefano M
OK, to jasne. Czy są jakieś wyraźne kryteria dla ? x1
Tim
1
@Tim nie ma wyraźnego progu cięcia, a odpowiedź zależy od precyzji implementacji zmiennoprzecinkowej (i rozwiązanego problemu). Chociaż expm1(x)powinna być dokładna do 1ULP w całym zakresie , stopniowo traci precyzję od kilku ULP, gdy x 1, do pełnego rozkładu, gdy x < ϵ , gdzie ϵ jest epsilonem maszynowym. 0x1exp(x) - 1x1x<ϵϵ
Stefano M
1

Aby rozwinąć różnicę między logi log1p, może być pomocne przywołanie wykresu, jeśli logarytm:

Logarithm

logx0ln(x)x0ln(x)ln(1e)=1ln(1e10)=10

x0ln(x+1)0ln(1+1e)0.31ln(1+1e10)0.000045log1p

1log01log1p

sfmiller940
źródło