Comparaison relative des nombres à virgule flottante

10

J'ai une fonction numérique f(x, y)renvoyant un double nombre à virgule flottante qui implémente une formule et je veux vérifier qu'elle est correcte par rapport aux expressions analytiques pour toutes les combinaisons de paramètres xet yque je suis intéressé. Quelle est la bonne façon de comparer les données calculées et nombres analytiques à virgule flottante?

Disons que les deux nombres sont aet b. Jusqu'à présent, je me suis assuré que les erreurs absolues ( abs(a-b) < eps) et relatives ( abs(a-b)/max(abs(a), abs(b)) < eps) sont inférieures à eps. De cette façon, il captera les inexactitudes numériques même si les chiffres sont disons autour de 1e-20.

Cependant, aujourd'hui, j'ai découvert un problème, la valeur numérique aet la valeur analytique bétaient:

In [47]: a                                                                     
Out[47]: 5.9781943146790832e-322

In [48]: b                                                                     
Out[48]: 6.0276008792632078e-322

In [50]: abs(a-b)                                                              
Out[50]: 4.9406564584124654e-324

In [52]: abs(a-b) / max(a, b)                                                  
Out[52]: 0.0081967213114754103

L'erreur absolue [50] est donc (évidemment) petite, mais l'erreur relative [52] est grande. J'ai donc pensé que j'avais un bug dans mon programme. En déboguant, je me suis rendu compte que ces nombres sont dénormaux . En tant que tel, j'ai écrit la routine suivante pour faire la comparaison relative appropriée:

real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
    r = 0
else
    r = d / m
end if
end function

Où tiny(1._dp)renvoie 2.22507385850720138E-308 sur mon ordinateur. Maintenant, tout fonctionne et j'obtiens simplement 0 comme erreur relative et tout va bien. En particulier, l'erreur relative ci-dessus [52] est erronée, elle est simplement causée par une précision insuffisante des nombres dénormaux. Ma mise en œuvre de la rel_errorfonction est-elle correcte? Dois-je simplement vérifier que abs(a-b)c'est moins que minuscule (= dénormal) et retourner 0? Ou devrais-je vérifier une autre combinaison, comme max(abs(a), abs(b))?

Je voudrais juste savoir quelle est la "bonne" façon.

floating-point Ondřej Čertík
la source

Réponses:

11

Vous pouvez vérifier directement les dénormals à l'aide isnormal()de math.h(C99 ou version ultérieure, POSIX.1 ou version ultérieure). Dans Fortran, si le module ieee_arithmeticest disponible, vous pouvez l'utiliser ieee_is_normal(). Pour être plus précis sur l'égalité floue, vous devez considérer la représentation en virgule flottante des dénormaux et décider de ce que vous voulez dire pour que les résultats soient suffisamment bons.

Plus précisément, pour croire que l'un ou l'autre résultat est correct, vous devez être sûr que vous n'avez pas perdu trop de chiffres à une étape intermédiaire. Le calcul avec des dénormals n'est généralement pas fiable et doit être évité en faisant redimensionner votre algorithme en interne. Pour vous assurer que votre mise à l'échelle interne a réussi, je recommande d'activer les exceptions à virgule flottante à l'aide feenableexcept()de C99 ou du ieee_arithmeticmodule dans Fortran.

Bien que votre application puisse capter le signal généré par les exceptions à virgule flottante, tous les noyaux que j'ai essayés ont réinitialisé l'indicateur matériel et fetestexcept()ne retournent donc pas de résultat utile. Lorsqu'ils sont exécutés avec -fp_trap, les programmes PETSc impriment (par défaut) une trace de pile lorsqu'une erreur en virgule flottante est déclenchée, mais n'identifient pas la ligne incriminée. Si vous exécutez dans un débogueur, le débogueur conserve l'indicateur matériel et interrompt l'expression incriminée. Vous pouvez vérifier la raison précise en appelant fetestexceptdepuis le débogueur où le résultat est un bit à bit OU des drapeaux suivants (les valeurs peuvent varier selon la machine, voir fenv.h; ces valeurs sont pour x86-64 avec glibc).

  • FE_INVALID = 0x1
  • FE_DIVBYZERO = 0x4
  • FE_OVERFLOW = 0x8
  • FE_UNDERFLOW = 0x10
  • FE_INEXACT = 0x20
Jed Brown
la source
Merci pour l'excellente réponse. L'expression analytique que je compare dans le régime asymptotique est exp(log_gamma(m+0.5_dp) - (m+0.5_dp)*log(t)) / 2pour m = 234, t = 2000. Il passe rapidement à zéro à mesure que j'augmente m. Tout ce que je veux m'assurer, c'est que ma routine numérique renvoie des nombres "corrects" (pour retourner zéro, c'est très bien aussi) à au moins 12 chiffres significatifs. Donc, si le calcul renvoie un nombre dénormal, il est simplement nul et il ne devrait pas y avoir de problème. Donc, juste la routine de comparaison doit être robuste contre cela.
Ondřej Čertík
5

Donald Knuth a une proposition pour un algorithme de comparaison à virgule flottante dans le volume 2 "Algorithmes séminariques" de "L'art de la programmation informatique". Il a été implémenté en C par Th. Belding (voir package fcmp ) et est disponible dans le GSL .

GertVdE
la source
2
Voici mon implémentation Fortran: gist.github.com/3776847 , notez que je dois de toute façon gérer explicitement les nombres dénormaux . Sinon, je pense que c'est à peu près l'équivalent de l'erreur relative, la seule différence est qu'au lieu de le faire abs(a-b)/max(a, b) < eps, nous le faisons abs(a-b)/2**exponent(max(a, b)) < eps, ce qui laisse pratiquement tomber la mantisse dans le max(a, b), donc à mon avis, la différence est négligeable.
Ondřej Čertík
5

Les nombres dénormalisés de manière optimale arrondie peuvent en effet avoir une erreur relative élevée. (Rincer cela à zéro tout en l'appelant une erreur relative est trompeur.)

Mais proche de zéro, le calcul d'erreurs relatives n'a aucun sens.

Par conséquent, avant même d'atteindre des nombres dénormalisés, vous devez probablement passer à la précision absolue (à savoir celle que vous souhaitez garantir dans ce cas).

yx|y−x|≤absacc+relacc∗max(|x|,|y|)

Les utilisateurs de votre code savent alors exactement combien de précision ils ont réellement.

Arnold Neumaier
la source
Êtes-vous sûr qu'il est inutile de calculer des erreurs relatives proches de zéro? Je pense que cela n'a de sens que s'il y a perte de précision (pour une raison quelconque). Si, par exemple, il y a une perte de précision pour x <1e-150 en raison de certains problèmes numériques (comme la soustraction de deux grands nombres), alors vous avez raison. Dans mon cas, cependant, les chiffres semblent exacts jusqu'à zéro, sauf lorsqu'ils atteignent les nombres dénormaux. Donc, dans mon cas, absacc = 1e-320 environ et je peux simplement vérifier abs(a-b) < tiny(1._dp)comme je le fais ci-dessus.
Ondřej Čertík
@ OndřejČertík: Dans ce cas, remplacez le 1e-150 par 1e-300 ou toute autre borne que vous pouvez vérifier. Dans tous les cas, très proche de zéro, vous faites une erreur absolue, et votre réclamation d'erreur devrait refléter cela plutôt que déclarer l'erreur relative comme nulle.
Arnold Neumaier,
Je vois. Je peux vérifier que tout fonctionne pour des nombres supérieurs à tiny(1._dp)=2.22507385850720138E-308(j'ai fait une erreur dans mon commentaire précédent, c'est 2e-308, pas 1e-320). C'est donc mon erreur absolue. Ensuite, je dois comparer l'erreur relative. Je vois votre point, je pense que vous avez raison. Merci!
Ondřej Čertík
1
|y−x|−absaccmax(|x|,|y|)