Petits résultats imprévisibles dans les exécutions d'un modèle déterministe

J'ai un modèle non négligeable (~ 5000 lignes) écrit en C. C'est un programme en série, sans génération de nombre aléatoire nulle part. Il utilise la bibliothèque FFTW pour les fonctions utilisant FFT - Je ne connais pas les détails de l'implémentation FFTW, mais je suppose que les fonctions y sont également déterministes (corrigez-moi si je me trompe).

Le problème que je ne peux pas comprendre est que j'obtiens de petites différences dans les résultats pour des exécutions identiques sur la même machine (même compilateur, mêmes bibliothèques).

J'utilise des variables de double précision, et pour sortir le résultat en variable valuepar exemple, j'émets: fprintf(outFID, "%.15e\n", value);ou
fwrite(&value, 1, sizeof(double), outFID);

Et j'obtiendrais constamment des différences telles que:
2.07843469652206 4 e-16 contre 2.07843469652206 3 e-16

J'ai passé beaucoup de temps à essayer de comprendre pourquoi. J'ai d'abord pensé qu'une de mes puces mémoire avait mal tourné, et je les ai commandées et remplacées, en vain. Par la suite, j'ai également essayé d'exécuter mon code sur la machine Linux d'un collègue, et j'obtiens des différences de même nature.

Qu'est-ce qui peut causer cela? C'est un petit problème maintenant, mais je me demande si c'est la "pointe de l'iceberg" (d'un problème grave).

Je pensais que je publierais ici au lieu de StackOverflow au cas où quelqu'un travaillant avec des modèles numériques aurait pu rencontrer ce problème. Si quelqu'un peut faire la lumière là-dessus, je serais très obligé.

Suivi des commentaires:
Christian Clason et Vikram: tout d'abord, merci de votre attention à ma question. Les articles que vous avez liés suggèrent que: 1. les erreurs d'arrondi limitent la précision, et 2. un code différent (tel que l'introduction d'instructions d'impression apparemment inoffensives) peut affecter les résultats jusqu'à la machine epsilon. Je dois préciser que je ne compare pas les effets fwriteet les fprintffonctions. J'utilise l'un OU l'autre. En particulier, le même exécutable est utilisé pour les deux exécutions. Je déclare simplement que le problème se produit si j'utilise fprintfOR fwrite.

Ainsi, le chemin du code (et l'exécutable) est le même et le matériel est le même. Avec tous ces facteurs externes maintenus constants, d'où vient l'aléatoire, fondamentalement? Je soupçonnais que le retournement de bits était dû à une mémoire défectueuse qui ne se retenait pas correctement, c'est pourquoi j'ai remplacé les puces de mémoire, mais cela ne semble pas être le problème ici, j'ai vérifié et vous l'avez indiqué. Mon programme génère des milliers de ces nombres en double précision en une seule fois, et il y a toujours une poignée aléatoire qui a des retournements de bits aléatoires.

Suivi du premier commentaire de Christian Clason: Pourquoi -il identique à 0 dans la précision de la machine? Le plus petit nombre positif pour un double est 2,22e-308, alors cela ne devrait-il pas être égal à 0? Mon programme génère des milliers de valeurs dans la gamme 10 ^ -16 (allant de 1e-15 à 8e-17) et nous avons vu des variations significatives dans notre projet de recherche, donc j'espère que nous n'avons pas regardé de manière absurde Nombres.$2\cdot 10^{-16}$

Suivi # 2 :
Ceci est un tracé de la série chronologique produite par le modèle, pour aider à la discussion de ramification dans les commentaires.

Il y a des aspects des systèmes informatiques modernes qui sont intrinsèquement non déterministes qui peuvent provoquer ce genre de différences. Tant que les différences sont très faibles par rapport à la précision requise de vos solutions, il n'y a probablement aucune raison de s'en inquiéter.

Un exemple de ce qui peut mal tourner d'après ma propre expérience. Considérons le problème du calcul du produit scalaire de deux vecteurs x et y.

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

$x_{i}y_{i}$

Par exemple, vous pouvez d'abord calculer le produit de deux vecteurs comme

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

puis comme

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

Comment cela pourrait-il arriver? Voici deux possibilités.

1. Calculs multithread sur coeurs parallèles. Les ordinateurs modernes ont généralement 2, 4, 8 ou même plus de cœurs de processeur qui peuvent fonctionner en parallèle. Si votre code utilise des threads parallèles pour calculer un produit scalaire sur plusieurs processeurs, toute perturbation aléatoire du système (par exemple, l'utilisateur a déplacé sa souris et l'un des cœurs de processeur doit traiter ce mouvement de la souris avant de revenir au produit scalaire) pourrait entraîner un changement dans l'ordre des ajouts.

2. ALignement des données et des instructions vectorielles. Les processeurs Intel modernes ont un ensemble spécial d'instructions qui peuvent fonctionner (par exemple) pour les nombres à virgule flottante à la fois. Ces instructions vectorielles fonctionnent mieux si les données sont alignées sur des limites de 16 octets. En règle générale, une boucle de produit scalaire diviserait les données en sections de 16 octets (4 flottants à la fois.) Si vous réexécutez le code une deuxième fois, les données pourraient être alignées différemment avec les blocs de mémoire de 16 octets afin que les ajouts soient effectué dans un ordre différent, résultant en une réponse différente.

Vous pouvez aborder le point 1 en faisant fonctionner votre code comme un seul thread et en désactivant tous les traitements parallèles. Vous pouvez aborder le point 2 en exigeant l'allocation de mémoire pour aligner les blocs de mémoire (en général, vous le feriez en compilant le code avec un commutateur tel que -align.) Si votre code donne toujours des résultats qui varient, il existe d'autres possibilités de recherche à.

Cette documentation d'Intel traite des problèmes pouvant entraîner une non reproductibilité des résultats avec la bibliothèque Intel Math Kernel. Un autre document d'Intel qui traite des commutateurs de compilateur à utiliser avec les compilateurs d'Intel.

La bibliothèque FFTW mentionnée peut s'exécuter en mode non déterministe.

Si vous utilisez le mode FFTW_MEASURE ou FFTW_PATIENT, les programmes vérifient au moment de l'exécution, quelles valeurs de paramètre fonctionnent le plus rapidement, puis utiliseront ces paramètres dans l'ensemble du programme. Parce que le temps d'exécution fluctuera évidemment un peu, les paramètres seront différents et le résultat des transformées de Fourier sera non déterministe. Si vous voulez une FFTW déterministe, utilisez le mode FFTW_ESTIMATE.

S'il est vrai que les changements d'ordre d'évaluation des termes d'expression peuvent très bien se produire en raison de scénarios de traitement multicœur / multithread, n'oubliez pas qu'il peut y avoir (même si c'est une longue perspective) une sorte de faille de conception matérielle au travail. Rappelez-vous le problème Pentium FDIV? (Voir https://en.wikipedia.org/wiki/Pentium_FDIV_bug ). Il y a quelque temps, j'ai travaillé sur un logiciel de simulation de circuits analogiques basé sur PC. Une partie de notre méthodologie consistait à développer des suites de tests de régression, que nous exécuterions contre les versions nocturnes du logiciel. Avec de nombreux modèles que nous avons développés, des méthodes itératives (par exemple Newton-Raphson ( https://en.wikipedia.org/wiki/Newton%27s_method)) et Runge-Kutta) ont été largement utilisés dans les algorithmes de simulation. Avec les appareils analogiques, il arrive souvent que des artefacts internes, tels que des tensions, des courants, etc., se produisent avec des valeurs numériques extrêmement petites. Ces valeurs, dans le cadre du processus de simulation, varient de façon incrémentielle dans le temps (simulé). L'amplitude de ces changements peut être extrêmement faible, et ce que nous avons souvent observé, c'est que les opérations FPU ultérieures sur de telles valeurs delta bordaient le seuil de "bruit" de la précision du FPU (le flottant 64 bits a une mantisse 53 bits, IIRC). Ceci, couplé au fait que nous devions souvent introduire du code de journalisation "PrintF" dans les modèles pour permettre le débogage (ah, les bons vieux jours!), Des résultats sporadiques pratiquement garantis, au quotidien! Et alors' s tout cela signifie? Vous devez vous attendre à voir des différences dans de telles circonstances, et la meilleure chose à faire est de définir et de mettre en œuvre un moyen de décider (ampleur, fréquence, tendance, etc.) quand / comment les ignorer.

Bien que l'arrondi à virgule flottante des opérations asynchrones puisse être le problème, je soupçonne que c'est quelque chose de plus banal. L'utilisation d'une variable non initialisée qui ajoute du caractère aléatoire à votre code par ailleurs déterministe. Il s'agit d'un problème courant qui est souvent ignoré par les développeurs car lorsque vous exécutez en mode débogage, toutes les variables sont initialisées à 0 lors de la déclaration. Lorsqu'elle n'est pas exécutée en mode débogage, la mémoire affectée à une variable a la valeur que la mémoire avait avant l'affectation. La mémoire n'est pas mise à zéro lors de l'affectation en tant qu'optimisation. Si cela se produit dans votre code, ce sera facile à corriger, moins dans le code des bibliothèques.