Względne porównanie liczb zmiennoprzecinkowych


10

Mam funkcję numeryczną f(x, y)zwracającą podwójną liczbę zmiennoprzecinkową, która implementuje pewną formułę i chcę sprawdzić, czy jest ona poprawna w stosunku do wyrażeń analitycznych dla wszystkich kombinacji parametrów xi yże jestem zainteresowany. Jaki jest właściwy sposób porównania obliczonego i analityczne liczby zmiennoprzecinkowe?

Powiedzmy, że dwie liczby to ai b. Do tej pory upewniłem się, że zarówno błędy bezwzględne ( abs(a-b) < eps), jak i względne ( abs(a-b)/max(abs(a), abs(b)) < eps) są mniejsze niż eps. W ten sposób wykryje niedokładności liczbowe, nawet jeśli liczby będą powiedzmy około 1e-20.

Jednak dzisiaj odkryłem problem, wartość liczbowa ai wartość analityczna bbyły następujące:

In [47]: a                                                                     
Out[47]: 5.9781943146790832e-322

In [48]: b                                                                     
Out[48]: 6.0276008792632078e-322

In [50]: abs(a-b)                                                              
Out[50]: 4.9406564584124654e-324

In [52]: abs(a-b) / max(a, b)                                                  
Out[52]: 0.0081967213114754103

Zatem błąd bezwzględny [50] jest (oczywiście) mały, ale błąd względny [52] jest duży. Pomyślałem więc, że mam błąd w moim programie. Podczas debugowania zdałem sobie sprawę, że liczby te są normalne . Jako taki napisałem następującą procedurę, aby wykonać właściwe porównanie względne:

real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
    r = 0
else
    r = d / m
end if
end function

Gdzie tiny(1._dp)zwraca 2.22507385850720138E-308 na moim komputerze. Teraz wszystko działa i po prostu dostaję 0 jako błąd względny i wszystko jest w porządku. W szczególności powyższy błąd względny [52] jest błędny, jest po prostu spowodowany niewystarczającą dokładnością liczb denormalnych. Czy moja implementacja rel_errorfunkcji jest poprawna? Czy powinienem po prostu sprawdzić, czy abs(a-b)jest mniejszy niż mały (= denormalny) i zwrócić 0? Czy powinienem sprawdzić inną kombinację, na przykład max(abs(a), abs(b))?

Chciałbym tylko wiedzieć, jaki jest „właściwy” sposób.

Odpowiedzi:


11

Można bezpośrednio sprawdzić denormals korzystających isnormal()z math.h(C99 lub nowszym POSIX.1 lub później). W module Fortran, jeśli moduł ieee_arithmeticjest dostępny, możesz użyć ieee_is_normal(). Aby być bardziej precyzyjnym na temat rozmytej równości, musisz wziąć pod uwagę zmiennoprzecinkową reprezentację denormali i zdecydować, co masz na myśli, aby wyniki były wystarczająco dobre.

Co więcej, aby uwierzyć, że którykolwiek z wyników jest poprawny, musisz mieć pewność, że nie straciłeś zbyt wielu cyfr na etapie pośrednim. Obliczenia z wartościami denormalnymi są generalnie zawodne i należy ich unikać poprzez wewnętrzną zmianę skali algorytmu. Aby upewnić się, że skalowanie wewnętrzne zakończyło się powodzeniem, zalecamy aktywację wyjątków zmiennoprzecinkowych za feenableexcept()pomocą C99 lub ieee_arithmeticmodułu w Fortran.

Mimo że aplikacja może przechwytywać sygnał podniesiony w wyjątkach zmiennoprzecinkowych, wszystkie jądra, które próbowałem zresetować flagę sprzętową, fetestexcept()nie zwracają użytecznego wyniku. Po uruchomieniu z -fp_trapprogramami PETSc (domyślnie) wydrukuje ślad stosu, gdy podniesiony zostanie błąd zmiennoprzecinkowy, ale nie zidentyfikuje linii obrażającej. Jeśli uruchomisz debugger, debuger zachowa flagę sprzętową i zepsuje obrażające wyrażenie. Możesz sprawdzić dokładny powód, dzwoniąc fetestexceptz debuggera, gdzie wynik jest bitowy LUB z następujących flag (wartości mogą się różnić w zależności od komputera, patrz fenv.h; te wartości dotyczą x86-64 z glibc).

  • FE_INVALID = 0x1
  • FE_DIVBYZERO = 0x4
  • FE_OVERFLOW = 0x8
  • FE_UNDERFLOW = 0x10
  • FE_INEXACT = 0x20

Dzięki za doskonałą odpowiedź. Wyrażenie analityczne, z którym porównuję w trybie asymptotycznym, jest exp(log_gamma(m+0.5_dp) - (m+0.5_dp)*log(t)) / 2dla m = 234, t = 2000. Wraz ze wzrostem szybko spada do zera m. Chcę tylko upewnić się, że moja procedura numeryczna zwraca „poprawne” liczby (powrót do zera też jest w porządku) do co najmniej 12 cyfr znaczących. Jeśli więc obliczenia zwracają liczbę normalną, to jest to po prostu zero i nie powinno być problemu. Dlatego właśnie procedura porównawcza musi być odporna na to.
Ondřej Čertík

5

Donald Knuth ma propozycję algorytmu porównywania zmiennoprzecinkowego w tomie 2 „Algorytmy seminumeryczne” z „The Art of Computer Programming”. Został zaimplementowany w C przez Th. Belding (patrz pakiet fcmp ) i jest dostępny w GSL .


2
Oto moja implementacja Fortrana : gist.github.com/3776847 , zauważ, że i tak muszę jawnie obsługiwać liczby normalne . W przeciwnym razie myślę, że to dość dużo odpowiednikiem błędu względnego, jedyną różnicą jest to, że zamiast robić abs(a-b)/max(a, b) < eps, robimy abs(a-b)/2**exponent(max(a, b)) < eps, co dość dużo po prostu spada mantysę W max(a, b), więc moim zdaniem różnica jest znikoma.
Ondřej Čertík

5

Optymalnie zaokrąglone liczby zdormalizowane mogą rzeczywiście mieć wysoki błąd względny. (Spłukanie go do zera, wciąż nazywając go błędem względnym, wprowadza w błąd).

Ale bliskie zeru, obliczanie względnych errosów jest bez znaczenia.

Dlatego nawet przed osiągnięciem liczb zdormalizowanych powinieneś prawdopodobnie przejść do absolutnej dokładności (mianowicie tej, którą chcesz zagwarantować w tym przypadku).

yx|yx|absacc+relaccmax(|x|,|y|)

Dzięki temu użytkownicy Twojego kodu dokładnie wiedzą, ile naprawdę mają dokładności.


Czy jesteś pewien, że obliczanie względnych błędów bliskich zeru nie ma sensu? Myślę, że nie ma to znaczenia, tylko jeśli nastąpi utrata dokładności (z jakiegokolwiek powodu). Jeśli na przykład występuje utrata dokładności dla x <1e-150 z powodu pewnych problemów liczbowych (takich jak odejmowanie dwóch dużych liczb), masz rację. W moim przypadku liczby wydają się jednak dokładne aż do zera, z wyjątkiem przypadków, gdy osiągają wartości normalne. Więc w moim przypadku absacc = 1e-320 lub więcej i mogę po prostu sprawdzić, abs(a-b) < tiny(1._dp)jak to robię powyżej.
Ondřej Čertík

@ OndřejČertík: W takim przypadku zastąp 1e-150 przez 1e-300 lub cokolwiek innego, co możesz zweryfikować. W każdym razie bardzo bliskie zeru popełniasz błąd bezwzględny, a twoja deklaracja błędu powinna to odzwierciedlać, a nie deklarować błąd względny jako zero.
Arnold Neumaier

Widzę. Mogę zweryfikować, że wszystkie działają dla liczb wyższych niż tiny(1._dp)=2.22507385850720138E-308(popełniłem błąd w poprzednim komentarzu, to 2e-308, a nie 1e-320). To jest mój absolutny błąd. Następnie muszę porównać błąd względny. Widzę twój punkt, myślę, że masz rację. Dzięki!
Ondřej Čertík,

1
@ OndřejČertík: Aby znaleźć dodatkowy błąd względny podany dla absacc, monitoruj maksimum . |yx|absaccmax(|x|,|y|)
Arnold Neumaier,
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.