Analiza błędu numerycznego w funkcji C ++

20

Załóżmy, że mam funkcję, która przyjmuje jako dane wejściowe kilka wartości zmiennoprzecinkowych (pojedyncze lub podwójne), wykonuje pewne obliczenia i generuje wyjściowe wartości zmiennoprzecinkowe (także pojedyncze lub podwójne). Pracuję przede wszystkim z MSVC 2008, ale planuję także współpracę z MinGW / GCC. Programuję w C ++.

Jaki jest typowy sposób programowego pomiaru ilości błędów w wynikach? Zakładając, że muszę użyć dowolnej biblioteki precyzji: jaka jest najlepsza taka biblioteka, jeśli nie zależy mi na szybkości?

użytkownik_123abc
źródło

Odpowiedzi:

17

Jeśli szukasz dobrego ograniczenia błędu zaokrąglania, niekoniecznie potrzebujesz biblioteki o precyzji aribtrary. Zamiast tego możesz użyć uruchomionej analizy błędów.

Nie udało mi się znaleźć dobrego odnośnika online, ale wszystko to opisano w rozdziale 3.3 książki Nicka Highama „Dokładność i stabilność algorytmów numerycznych”. Pomysł jest dość prosty:

  1. Ponownie wprowadź współczynnik do kodu, aby w każdym wierszu było jedno przypisanie jednej operacji arytmetycznej.
  2. Dla każdej zmiennej np. xUtwórz zmienną, x_errktóra jest inicjalizowana na zero, gdy xzostanie przypisana stała.
  3. Dla każdej operacji, np. z = x * yZaktualizuj zmienną z_errprzy użyciu standardowego modelu arytmetyki zmiennoprzecinkowej oraz wynikających zz niej błędów błędów x_erri y_err.
  4. Zwracana wartość twojej funkcji powinna wtedy również mieć przypisaną odpowiednią _errwartość. Jest to zależne od danych ograniczenie całkowitego błędu zaokrąglenia.

Trudna część to krok 3. W przypadku najprostszych operacji arytmetycznych można zastosować następujące reguły:

  • z = x + y -> z_err = u*abs(z) + x_err + y_err
  • z = x - y -> z_err = u*abs(z) + x_err + y_err
  • z = x * y -> z_err = u*abs(z) + x_err*abs(y) + y_err*abs(x)
  • z = x / y -> z_err = u*abs(z) + (x_err*abs(y) + y_err*abs(x))/y^2
  • z = sqrt(x) -> z_err = u*abs(z) + x_err/(2*abs(z))

gdzie u = eps/2jest zaokrąglenie jednostki. Tak, zasady +i -są takie same. Reguły dla każdej innej operacji op(x)można łatwo wyodrębnić za pomocą rozszerzenia serii Taylora zastosowanego wyniku op(x + x_err). Lub możesz spróbować google. Lub korzystając z książki Nicka Highama.

Jako przykład rozważmy następujący kod Matlab / Octave, który ocenia wielomiany we współczynnikach aw punkcie xprzy użyciu schematu Hornera:

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        s = a(k) + x*s;
    end

W pierwszym etapie podzieliliśmy dwie operacje na s = a(k) + x*s:

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        z = x*s;
        s = a(k) + z;
    end

Następnie wprowadzamy _errzmienne. Zauważ, że dane wejściowe ai xzakłada się, że są dokładne, ale równie dobrze możemy również wymagać od użytkownika przekazania odpowiednich wartości dla a_erri x_err:

function [ s , s_err ] = horner ( a , x )
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = ...;
        s = a(k) + z;
        s_err = ...;
    end

Na koniec stosujemy zasady opisane powyżej, aby uzyskać warunki błędów:

function [ s , s_err ] = horner ( a , x )
    u = eps/2;
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = u*abs(z) + s_err*abs(x);
        s = a(k) + z;
        s_err = u*abs(s) + z_err;
    end

Zauważ, że ponieważ nie mamy a_errlub x_errnp. Przyjmuje się, że są zerowe, odpowiednie wyrażenia są po prostu ignorowane w wyrażeniach błędów.

Zrobione! Mamy teraz schemat Hornera, który zwraca wynik szacunku błędu zależnego od danych (uwaga: jest to górna granica błędu) obok wyniku.

Na marginesie, skoro używasz C ++, możesz rozważyć utworzenie własnej klasy dla wartości zmiennoprzecinkowych, która przenosi ten _errtermin i przeciąża wszystkie operacje arytmetyczne, aby zaktualizować te wartości, jak opisano powyżej. W przypadku dużych kodów może to być łatwiejsza, choć mniej obliczeniowa, trasa. Powiedziawszy to, możesz być w stanie znaleźć taką klasę online. Szybkie wyszukiwanie w Google dało mi ten link .

±ux(1±u)

Pedro
źródło
1
+1 za tę analizę, ponieważ jest interesująca. Lubię pracę Highama. Niepokoi mnie to, że wymaganie od użytkownika ręcznego napisania tego dodatkowego kodu (zamiast półautomatycznych, takich jak arytmetyka interwałowa) może być podatne na błędy, ponieważ liczba operacji numerycznych staje się duża.
Geoff Oxberry
1
@GeoffOxberry: Całkowicie zgadzam się z problemem złożoności. W przypadku większych kodów zdecydowanie zaleciłbym napisanie klasy / typu danych, który przeciąża operacje na podwójnych kartach, tak aby każda operacja była poprawnie wykonywana tylko raz. Jestem dość zaskoczony, że wydaje się, że nie ma czegoś takiego dla Matlab / Octave.
Pedro
Podoba mi się ta analiza, ale skoro obliczenia warunków błędu są również wykonywane w liczbach zmiennoprzecinkowych, czy te warunki błędów nie będą nieprecyzyjne z powodu błędów zmiennoprzecinkowych?
plasmacel,
8

Fajną przenośną i otwartą biblioteką dla dowolnej precyzji arytmetyki zmiennoprzecinkowej (i wielu innych) jest NTL Victora Shoupa , która jest dostępna w formie źródła C ++.

Na niższym poziomie znajduje się biblioteka Bignum GNU Multiple Precision (GMP) , również pakiet open source.

NTL może być używany z GMP, ponieważ wymagana jest szybsza wydajność, ale NTL zapewnia własne podstawowe procedury, które z pewnością są użyteczne, jeśli „nie zależy ci na szybkości”. GMP twierdzi, że jest „najszybszą biblioteką bignum”. GMP jest w dużej mierze napisany w C, ale ma interfejs C ++.

Dodano: Chociaż arytmetyka interwałów może określać górną i dolną granicę dokładnej odpowiedzi w sposób zautomatyzowany, nie mierzy ona dokładnie błędu w „standardowym” precyzyjnym obliczeniu, ponieważ rozmiar interwału zwykle rośnie z każdą operacją (względną lub błąd bezwzględny).

Typowym sposobem znalezienia wielkości błędu, zarówno w przypadku błędów zaokrąglania, jak i błędów dyskretyzacji itp., Jest obliczenie dodatkowej wartości dokładności i porównanie jej ze „standardową” wartością dokładności. Tylko niewielka dodatkowa precyzja jest potrzebna do określenia samego rozmiaru błędu do rozsądnej dokładności, ponieważ same błędy zaokrąglania są znacznie większe w „standardowej” precyzji niż w obliczeniach o dodatkowej precyzji.

Można to zilustrować poprzez porównanie obliczeń o pojedynczej i podwójnej precyzji. Zauważ, że w C ++ wyrażenia pośrednie są zawsze obliczane z (przynajmniej) podwójną precyzją, więc jeśli chcemy zilustrować, jak wyglądałoby obliczenie w „czystej” pojedynczej precyzji, musimy przechowywać wartości pośrednie w pojedynczej precyzji.

Fragment kodu C.

    float fa,fb;
    double da,db,err;
    fa = 4.0;
    fb = 3.0;
    fa = fa/fb;
    fa -= 1.0;

    da = 4.0;
    db = 3.0;
    da = da/db;
    da -= 1.0;

    err = fa - da;
    printf("Single precision error wrt double precision value\n");
    printf("Error in getting 1/3rd is %e\n",err);
    return 0;

Dane wyjściowe z góry (łańcuch narzędzi Cygwin / MinGW32 GCC):

Single precision error wrt double precision value
Error in getting 1/3rd is 3.973643e-08

Zatem błąd dotyczy tego, czego oczekuje się w zaokrąglaniu 1/3 do pojedynczej precyzji. One nie (ja podejrzewam) dbałość o coraz więcej niż kilka miejsc po przecinku w błędzie poprawne, ponieważ pomiar błędu jest do wielkości, a nie dla ścisłości.

matematyka
źródło
Twoje podejście jest zdecydowanie matematyczne. Myślę, że kompromis jest rygorystyczny; ludzie, którzy są pedantyczni w kwestii błędu, wskażą na rygor arytmetyki interwałów, ale podejrzewam, że w wielu aplikacjach obliczenia z dodatkową precyzją byłyby wystarczające, a wynikowe oszacowania błędów prawdopodobnie byłyby dokładniejsze, jak wskazałeś.
Geoff Oxberry
Właśnie takie podejście sobie wyobrażałem. Mógłbym wypróbować kilka z tych różnych technik, aby zobaczyć, które są najbardziej odpowiednie dla mojej aplikacji. Aktualizacja przykładowego kodu jest bardzo mile widziana!
user_123abc
7

GMP (tj. Biblioteka GNU Multiple Precision) jest najlepszą biblioteką o dowolnej precyzji, jaką znam.

Nie znam żadnego programowego sposobu pomiaru błędu w wynikach dowolnej funkcji zmiennoprzecinkowej. Jedną z rzeczy, które możesz wypróbować, jest obliczenie przedłużenia interwału funkcji za pomocą arytmetyki interwałów . W C ++ trzeba by użyć biblioteki do obliczania rozszerzeń przedziałów; jedną z takich bibliotek jest biblioteka arytmetyczna przedziału czasu doładowania. Zasadniczo, aby zmierzyć błąd, podajemy jako argumenty interwały funkcji, które mają szerokość 2-krotności jednostkowego zaokrąglenia (z grubsza), wyśrodkowaną na interesujących wartościach, a następnie wynikiem byłby zbiór przedziałów, szerokość co dałoby pewne zachowawcze oszacowanie błędu. Trudność z tym podejściem polega na tym, że arytmetyka przedziałowa stosowana w ten sposób może przeceniać błąd o znaczne kwoty, ale takie podejście jest najbardziej „programowe”, jakie mogę sobie wyobrazić.

Geoff Oxberry
źródło
Ach, właśnie zauważyłem arytmetykę przedziałową wspomnianą w twojej odpowiedzi ... Głosowałem!
Ali
2
Richard Harris napisał doskonałą serię artykułów w dzienniku ACCU Overload o Floating Point Blues . Jego artykuł na temat arytmetyki przedziałowej znajduje się w Przeciążeniu 103 ( pdf , p19-24).
Mark Booth,
6

Rygorystyczne i automatyczne oszacowanie błędu można osiągnąć za pomocą analizy interwałowej . Pracujesz w odstępach zamiast liczb. Na przykład dodanie:

[a,b] + [c,d] = [min(a+c, a+d, b+c, b+d), max (a+c, a+d, b+c, b+d)] = [a+c, b+d]

Zaokrąglanie można również przeprowadzić rygorystycznie, patrz Arytmetyka zaokrąglonych interwałów .

O ile dane wejściowe składają się z wąskich przedziałów, szacunki są OK i są naprawdę tanie do obliczenia. Niestety błąd jest często przeceniany, patrz problem zależności .

Nie znam żadnej biblioteki arytmetycznej o dowolnym przedziale precyzji.

To, czy arytmetyka interwałowa może zaspokoić twoje potrzeby, zależy od twojego problemu.

Ali
źródło
4

Biblioteka GNU MPFR jest biblioteką pływak arbitralny precyzji, który ma wysoką dokładność (w szczególności prawidłowego zaokrąglenia dla wszystkich operacji, które nie jest tak proste, jak to brzmi) jako jeden z ich głównych punktów ostrości. Używa GNU MP pod maską. Ma rozszerzenie o nazwie MPFI, które wykonuje arytmetykę interwałów, co - jak sugeruje odpowiedź Geoffa - może się przydać do celów weryfikacji: zwiększaj precyzję działania, aż wynikowy interwał znajdzie się w niewielkich granicach.

Jednak to nie zawsze zadziała; w szczególności niekoniecznie jest to skuteczne, jeśli wykonujesz coś takiego jak integracja numeryczna, gdzie każdy krok niesie „błąd” niezależny od problemów z zaokrąglaniem. W takim przypadku wypróbuj specjalistyczny pakiet, taki jak COSY infinity, który robi to bardzo dobrze, używając specyficznych algorytmów w celu ograniczenia błędu integracji (i stosując tak zwane modele Taylora zamiast przedziałów).

Erik P.
źródło
Zgadzam się; integracja numeryczna to zdecydowanie przypadek, w którym arytmetyka przedziałów naiwnych może pójść nie tak. Jednak nawet modele Taylora używają arytmetyki przedziałowej. Znam pracę Makino i Berza i wierzę, że używają modelu Taylora w sensie RE Moore, chociaż używają także sztuczek, które nazywają „algebrą różnicową”.
Geoff Oxberry
@GeoffOxberry: Tak - myślę, że ta algebra różniczkowa jest wystarczająca do ograniczenia błędu na etapie integracji.
Erik P.
0

Powiedziano mi, że MPIR jest dobrą biblioteką do użycia, jeśli pracujesz z Visual Studio:

http://mpir.org/

rybak
źródło
Witamy w SciComp.SE! Czy mógłbyś dodać kilka szczegółów na temat tego, jak można wykorzystać tę bibliotekę do pomiaru błędu obliczeń zmiennoprzecinkowych?
Christian Clason
Spróbuję; Właściwie nie skonfigurowałem jeszcze MPIR na moim komputerze! Skonfigurowałem GMP i MPFR.
fishermanhat