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ę to zajęło, ale udało mi się dowiedzieć, z czego wynika 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?