Dokładność numeryczna

Poniższy wpis autorstwa pana Grzegorza Knora otrzymał wyróżnienie w naszym czerwcowym konkursie

Wstęp

Poniższy wpis bazuje na wątku dokładność numeryczna z polskiego forum matlab.pl.
Otóż w 2009 roku zetknąłem sie z problemem, który może być zilustrowany za pomocą następującego pseudokodu:

Dla każdego k ze zbioru [0,1, 0,2, 0,3,... 5]

Jeśli zaokrąglona wartość k jest równa k

Wydrukuj k


Czyli chodzi o wybranie liczb całkowitych z naszego zbioru. Rozwiązanie tego zadania w MATLABie jest bardzo proste:

k = 0.1:0.1:5;% wektor k [0,1, 0,2, 0,3,... 5]
idx = k == round(k); % szukanie wartości całkowitych
disp(k(idx)) % drukowanie wyniku

1 2 3 4 5


Jak widać rozwiązanie jest oczywiste, to liczby 1, 2, 3, 4 i 5. Ale rozbudujmy nieco nasze zadanie i załóżmy, że musimy wykonać pewne obliczenia w pętli z krokiem co 0,1 sekundy, zaś wyniki (np. wykres) chcemy drukować dla pełnych sekund. Ponownie, rozwiązanie nie jest skomplikowane:

for k = 0.1:0.1:5
    % skomplikowane obliczenia
    if k == round(k)
        % drukuj wynik dla całkowitych wartości k
        disp(k)
    end
end

1
2
4
5

Wynik jest jednak zaskakujący. Co się stało z wynikami w 3 sekundzie? Zbadajmy to poprzez małą modyfikację powyższego kodu. Liczby całkowite w naszym wektorze są na pozycjach 10, 20, 30, 40 i 50. Wydrukujmy zatem różnice pomiędzy liczbami, które powinny być całkowite i zaokrąglonymi wartościami:

idx = 0;
for k = 0.1:0.1:5
    idx = idx + 1;
    if mod(idx,10)==0  
        fprintf('liczba %i, roznica: %E',round(k),k-round(k))
    end
end

liczba 1, roznica: 0.000000E+00
liczba 2, roznica: 0.000000E+00
liczba 3, roznica: 4.440892E-16
liczba 4, roznica: 0.000000E+00
liczba 5, roznica: 0.000000E+00


W istocie, w przypadku liczby 3 różnica nie jest zerowa. Dlaczego? Jest to problem z reprezentacją liczb zmiennoprzecinkowych, a konkretnie liczba 0,1 (którą MATLAB dodaje do wartości k w każdym kroku pętli) nie może być zapisana w postaci binarnej (więcej informacji można znaleźć w Wikipedii lub dokumentacji MATLABa).

Prawdziwy problem

 

Wytłumaczenie powyższego problemu reprezentacją liczb zmiennoprzecinkowych jest zadowalające. Ale spróbujmy nieco zmodyfikować wyjściowy kod i stwórzmy wektor k przed pętla. Jeśli pętla będzie przebiegała po wszystkich wartościach k wynik powinien pozostać taki sam:

k= 0.1:0.1:5;
for m=1:length(k)
    % skomplikowane obliczenia
    if k(m) == round(k(m))
        % drukuj wynik dla całkowitych wartości k
        disp(k(m))
    end
end

Tym razem wynik jest poprawny, ale nie tego oczekiwaliśmy. Na szczęście, wsparcie techniczne MATLABa działa całkiem sprawnie 🙂 Poniżej odpowiedź administratora forum mpi wyjaśniająca różnice (przy czym w oryginalnym pytaniu k oznaczało wektor tworzony przed pętlą, zaś k1 było iteratorem pętli; w tym wpisie zarówno wektor tworzony przed pętlą jak i iterator oznaczane są wspólnym symbolem k):
Trochę mi to zajęło, ale udało mi się dowiedzieć skąd się bierze różnica. Wektor k jest tworzony i liczony przez BLAS (część jądra obliczeniowego Matlaba), który trzyma ją w 80 bitowych specjalnych rejestrach, żeby było dokładniej. Wektor k1 jest tylko iteratorem pętli, więc ze względu na oszczędność pamięci jest tworzony i liczony w zwykłej zmiennej więc jest mniej dokładny. Pamiętajcie, że producent NIE zaleca stosowania pętli for ze względu na tzw. małą elastyczność. Taką pętlę bardzo ciężko zoptymalizować, w przeciwieństwie do operacji wektorowych. Z operacjami wektorowymi JIT potrafi wyprawiać prawdziwe cuda. Pętlę for uważa się za zło konieczne, a już nigdy, przenigdy nie powinno się wkładać precyzyjnych operacji matematycznych do iteratora.
Jako komentarz do tej odpowiedzi dodam tylko, że od 2009 roku trochę się zmieniło i pętle for nie są już takie wolne w MATLABie.

Diabeł tkwi w szczegółach

A jak się dowiedzieć kiedy stworzony wektor będzie używał BLAS, a kiedy nie? Przyznam, że nie znam odpowiedzi na to pytanie, i nie potrafię jej znaleźć w dokumentacji MATLABa. Co więcej czasem MATLAB mnie zaskakuje. Dla przykładu raz jeszcze zmodyfikujmy wyjściową pętlę. Tym razem zmiana będzie kosmetyczna, dodajmy nawiasy kwadratowe. A jaki będzie wynik?

for k = [0.1:0.1:5]
    % skomplikowane obliczenia
    if k == round(k)
        % drukuj wynik dla calkowitych wartosci k
        disp(k)
    end
end

1
2
3
4
5

Dodanie nawiasów kwadratowych powoduje, że tym razem iterator pętli jest bardziej dokładny niż w przypadku bez nawiasów, mimo iż MATLAB (a konkretnie Code Analyzer) doradza nam, że użycie nawiasów jest zbędne w tym przypadku:

Powyższy przykład pokazuje, iż żartobliwe powiedzenie, że programowanie jest dziedziną eksperymentalną ma w sobie ziarno prawdy 🙂

Deser

Czy są jeszcze inne przykłady podobnego typu? W swojej pracy zetknąłem się z podobnym zagadnieniem podczas debugowania kodu w MATLABie 2010a. Mianowicie wynik obliczeń różnił się w zależności od tego czy funkcja była wywołana zwyczajnie, czy w trybie debuggera (z kropką do zatrzymania się w wybranej linijce). Jeśli stworzymy przykładową funkcje foo:

function out = foo(a,b)
out = a < b;
end

To porównanie liczb bardzo bliskich siebie może dawać różny wynik przy zwykłym uruchomieniu i inny podczas debugowania. Przyczyna była podobna jak w przypadku iteratora pętli, MATLAB wykonywał obliczenia z inną dokładnością w trybie debuggera. Na szczęście problem został juz naprawiony i nie występuje w obecnej wersji MATLABa.
A czy Wy natknęliście się na problemy podobnego typu w swojej pracy z MATLABem?

2 myśli nt. „Dokładność numeryczna”

  1. Bardzo ciekawie rozpisana sprawa. Ja potrzebowałem chyba godziny, aby znaleźć taki mankament u siebie:

    Poszukiwałem najlepszej wartości pewnego parametru z zakresu od 0 do 1 z krokiem 0.1 i oczywiście robiłem to w pętli for ,,z dwukropkiem''.

    Program odpowiada np. że dla wartości parametru równej 0.3 znalazł najlepszy wynik. Została ona tak wypisana na ekran.

    I kiedy potem ręcznie ustawiam wartość parametru na 0.3 wpisując ją z klawiatury, nie otrzymuję tych samych wyników.

    Poniższy kod opisuje mój problem w najbardziej trywialny sposób:

    i = 0:0.1:1;
    j = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1];
    disp(i==j)

    W jego wyniku otrzymałem:
    1 1 1 0 1 1 1 1 1 1 1

    Podpowiem, że testowałem ten problem też na innych ,,konkurencyjnych'' środowiskach i było jeszcze więcej zer, szczególnie przy definiowaniu innych zakresów i innych kroków w pętli.
    Pozdrawiam

Dodaj komentarz

Twój adres email nie zostanie opublikowany. Pola, których wypełnienie jest wymagane, są oznaczone symbolem *