Małe, nieprzewidywalne wyniki w seriach modelu deterministycznego

Mam spory model (~ 5000 linii) napisany w C. Jest to program szeregowy, bez generowania liczb losowych w dowolnym miejscu. Wykorzystuje bibliotekę FFTW dla funkcji korzystających z FFT - nie znam szczegółów implementacji FFTW, ale zakładam, że funkcje w niej również są deterministyczne (poprawcie mnie, jeśli się mylę).

Problem, którego nie rozumiem, polega na tym, że otrzymuję małe różnice w wynikach dla identycznych uruchomień na tym samym komputerze (ten sam kompilator, te same biblioteki).

Używam zmiennych o podwójnej precyzji i na przykład w wyniku wypisuję zmienną value: fprintf(outFID, "%.15e\n", value);lub
fwrite(&value, 1, sizeof(double), outFID);

I ciągle otrzymywałbym różnice, takie jak:
2.07843469652206 4 e-16 vs. 2.07843469652206 3 e-16

Spędziłem dużo czasu próbując zrozumieć, dlaczego tak jest. Początkowo myślałem, że jeden z moich układów pamięci zepsuł się i zamówiłem je i wymieniłem, ale bezskutecznie. Następnie spróbowałem uruchomić mój kod na maszynie Linuksa kolegi i dostaję różnice o tym samym charakterze.

Co może być tego przyczyną? Jest to teraz niewielka kwestia, ale zastanawiam się, czy jest to „wierzchołek góry lodowej” (poważny problem).

Myślałem, że opublikuję tutaj zamiast StackOverflow na wypadek, gdyby ktoś pracujący z modelami numerycznymi napotkał ten problem. Gdyby ktoś mógł rzucić na to światło, byłbym bardzo zobowiązany.

Kontynuacja komentarzy:
Christian Clason i Vikram: po pierwsze, dziękuję za uwagę na moje pytanie. Artykuły, które podłączyłeś, sugerują, że: 1. błędy zaokrąglania ograniczają dokładność, oraz 2. inny kod (np. Wprowadzanie pozornie nieszkodliwych instrukcji drukowania) może wpływać na wyniki aż do epsilon maszyny. Powinienem wyjaśnić, że nie porównuję efektów fwritei fprintffunkcji. Używam jednego LUB drugiego. W szczególności ten sam plik wykonywalny jest używany dla obu przebiegów. Po prostu stwierdzam, że problem występuje bez względu na to, czy używam fprintfOR fwrite.

Zatem ścieżka do kodu (i wykonywalna) jest taka sama, a sprzęt jest taki sam. Przy tych wszystkich czynnikach zewnętrznych utrzymywanych na stałym poziomie, skąd zasadniczo bierze się przypadkowość? Podejrzewałem, że bit się odwrócił, ponieważ uszkodzona pamięć nie zachowuje się trochę poprawnie, dlatego wymieniłem układy pamięci, ale to nie wydaje się być problemem tutaj, zweryfikowałem i wskazałeś. Mój program generuje tysiące tych podwójnie precyzyjnych liczb w jednym przebiegu i zawsze istnieje losowa garść, która ma losowe odwrócenie bitów.

Kontynuacja pierwszego komentarza Christiana Clasona: Dlaczego tym samym co 0 w zakresie precyzji maszyny? Najmniejsza liczba dodatnia dla podwójnego to 2.22e-308, więc czy nie powinno to być równe 0? Mój program generuje tysiące wartości w zakresie 10 ^ -16 (od 1e-15 do 8e-17) i widzieliśmy znaczące różnice w naszym projekcie badawczym, więc mam nadzieję, że nie patrzyliśmy na nonsensowne liczby.$2\cdot 10^{-16}$

Kontynuacja # 2 :
Jest to wykres szeregów czasowych generowanych przez model, aby pomóc w dyskusji na temat odejścia w komentarzach.

Istnieją aspekty współczesnych systemów obliczeniowych, które z natury są niedeterministyczne, które mogą powodować tego rodzaju różnice. Dopóki różnice są bardzo małe w porównaniu z wymaganą dokładnością twoich rozwiązań, prawdopodobnie nie ma powodu, aby się tym martwić.

Przykład tego, co może pójść nie tak, na podstawie własnego doświadczenia. Rozważ problem obliczania iloczynu kropkowego dwóch wektorów x i y.

$d=\sum_{i=1}^{n} x_{i}y_{i}$

Obliczenie tego produktu kropkowego wymaga obliczenia produktów a następnie zsumowania wyników. Mnożenie zmiennoprzecinkowe powinno za każdym razem dawać dokładnie takie same wyniki. Jeśli uzupełnienia zmiennoprzecinkowe są obliczane za każdym razem w tej samej kolejności, suma powinna być identyczna. Ponieważ jednak dodawanie zmiennoprzecinkowe nie jest asocjacyjne, możesz uzyskać różne wyniki, jeśli iloczyn dwóch wektorów jest obliczany w taki sposób, że dodawanie jest wykonywane w innej kolejności. $x_{i}y_{i}$

Na przykład możesz najpierw obliczyć iloczyn dwóch wektorów jako

$d=((x_{1}y_{1})+(x_{2}y_{2}))+(x_{3}y_{3})$

a następnie jako

$d=(x_{1}y_{1})+((x_{2}y_{2})+(x_{3}y_{3}))$ .

Jak to się mogło stać? Oto dwie możliwości.

1. Obliczenia wielowątkowe na rdzeniach równoległych. Nowoczesne komputery zwykle mają 2, 4, 8 lub nawet więcej rdzeni procesorów, które mogą pracować równolegle. Jeśli kod używa równoległych wątków do obliczania iloczynu na wielu procesorach, wówczas dowolne zaburzenie systemu (np. Użytkownik przesunął mysz i jeden z rdzeni procesora musi przetworzyć ten ruch myszy przed powrotem do iloczynu) spowodować zmianę kolejności dodawania.

2. Wyrównanie danych i instrukcji wektorowych. Nowoczesne procesory Intel mają specjalny zestaw instrukcji, które mogą działać (na przykład) dla liczb zmiennoprzecinkowych jednocześnie. Te instrukcje wektorowe działają najlepiej, jeśli dane są wyrównane do 16 bajtów granic. Zwykle pętla iloczynu dzieli dane na sekcje po 16 bajtów (4 zmienne na raz.) Jeśli ponownie uruchomisz kod po raz drugi, dane mogą zostać wyrównane w inny sposób z 16-bajtowymi blokami pamięci, dzięki czemu dodawane są wykonywane w innej kolejności, co skutkuje inną odpowiedzią.

Możesz rozwiązać punkt 1, uruchamiając kod jako pojedynczy wątek i wyłączając wszystkie przetwarzanie równoległe. Możesz zająć się punktem 2, wymagając alokacji pamięci, aby wyrównać bloki pamięci (zwykle robisz to, kompilując kod za pomocą przełącznika, takiego jak -align.) Jeśli kod nadal daje różne wyniki, istnieją inne możliwości wyszukiwania w.

Ta dokumentacja firmy Intel omawia problemy, które mogą prowadzić do braku powtarzalności wyników w bibliotece Intel Math Kernel Library. Kolejny dokument firmy Intel omawiający przełączniki kompilatora do użycia z kompilatorami Intela.

Wspomniana biblioteka FFTW może działać w trybie niedeterministycznym.

Jeśli używasz trybu FFTW_MEASURE lub FFTW_PATIENT, programy sprawdzają w czasie wykonywania, które wartości parametrów działają najszybciej, a następnie będą używać tych parametrów w całym programie. Ponieważ czas pracy oczywiście nieco się waha, parametry będą różne, a wynik przekształceń Fouriera będzie niedeterministyczny. Jeśli chcesz deterministyczne FFTW, użyj trybu FFTW_ESTIMATE.

Chociaż prawdą jest, że zmiany kolejności oceny wyrażeń mogą bardzo dobrze wystąpić ze względu na scenariusze przetwarzania wielordzeniowego / wielowątkowego, nie zapominaj, że może istnieć (nawet jeśli jest to długa szansa) jakaś usterka w projektowaniu sprzętu. Pamiętasz problem Pentium FDIV? (Zobacz https://en.wikipedia.org/wiki/Pentium_FDIV_bug ). Jakiś czas temu pracowałem nad oprogramowaniem komputerowym do symulacji obwodów analogowych. Część naszej metodologii obejmowała opracowanie pakietów testów regresji, które działałyby przeciwko nocnym wersjom oprogramowania. W przypadku wielu opracowanych modeli metody iteracyjne (np. Newton-Raphson ( https://en.wikipedia.org/wiki/Newton%27s_method) i Runge-Kutta) były szeroko stosowane w algorytmach symulacyjnych. W przypadku urządzeń analogowych często zdarza się, że artefakty wewnętrzne, takie jak napięcia, prądy itp., Mają wyjątkowo małe wartości liczbowe. Wartości te, jako część procesu symulacji, zmieniają się stopniowo w czasie (symulowanym). Skala tych zmian może być bardzo mała i często obserwowaliśmy, że kolejne operacje FPU na takich wartościach delta graniczyły z progiem „szumu” precyzji FPU (pływanie 64-bitowe ma 53-bitową mantysę, IIRC). To, w połączeniu z faktem, że często musieliśmy wprowadzać do modeli kod rejestrujący „PrintF”, aby umożliwić debugowanie (ach, dobre stare dni!), Praktycznie gwarantowało sporadyczne wyniki, na co dzień! Więc co' czy to wszystko znaczy? W takich okolicznościach należy spodziewać się różnic, a najlepszą rzeczą jest zdefiniowanie i wdrożenie sposobu decydowania (wielkość, częstotliwość, trend itp.), Kiedy / jak je zignorować.

Podczas gdy zaokrąglanie zmiennoprzecinkowe z operacji asynchronicznych może być problemem, podejrzewam, że jest to coś bardziej banalnego. Zastosowanie niezainicjowanej zmiennej, która dodaje losowość do twojego kodu deterministycznego. Jest to częsty problem, który jest często pomijany przez programistów, ponieważ podczas uruchamiania w trybie debugowania wszystkie zmienne są inicjowane na 0 w momencie deklaracji. Gdy pamięć nie jest uruchomiona w trybie debugowania, pamięć przypisana do zmiennej ma dowolną wartość, jaką pamięć miała przed przypisaniem. Pamięć nie jest zerowana przy przypisaniu jako optymalizacja. Jeśli dzieje się tak w twoim kodzie, łatwo będzie to naprawić, w mniejszym stopniu w kodzie bibliotek.