Analiza błędu numerycznego w funkcji C ++

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?

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 .

$\pm u$$x(1 \pm u)$

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.

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ć.

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.

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).

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

http://mpir.org/