Relativer Vergleich von Gleitkommazahlen

10

Ich habe eine numerische Funktion, f(x, y)die eine doppelte Gleitkommazahl zurückgibt, die eine Formel implementiert, und ich möchte überprüfen, ob sie gegen analytische Ausdrücke für alle Kombinationen der Parameter korrekt ist xund ymich interessiert. Was ist der richtige Weg, um die berechneten und zu vergleichen? analytische Gleitkommazahlen?

Nehmen wir an, die beiden Zahlen sind aund b. Bisher habe ich sichergestellt, dass sowohl absolute ( abs(a-b) < eps) als auch relative ( abs(a-b)/max(abs(a), abs(b)) < eps) Fehler kleiner als eps sind. Auf diese Weise werden numerische Ungenauigkeiten erfasst, selbst wenn die Zahlen etwa 1e-20 betragen.

Heute habe ich jedoch ein Problem entdeckt, der numerische Wert aund der analytische Wert bwaren:

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

Der absolute Fehler [50] ist also (offensichtlich) klein, aber der relative Fehler [52] ist groß. Also dachte ich, dass ich einen Fehler in meinem Programm habe. Durch das Debuggen wurde mir klar, dass diese Zahlen denormal sind . Als solches habe ich die folgende Routine geschrieben, um den richtigen relativen Vergleich durchzuführen:

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

Wo tiny(1._dp)gibt 2.22507385850720138E-308 auf meinem Computer zurück. Jetzt funktioniert alles und ich bekomme einfach 0 als relativen Fehler und alles ist in Ordnung. Insbesondere ist der obige relative Fehler [52] falsch, er wird einfach durch eine unzureichende Genauigkeit der denormalen Zahlen verursacht. Ist meine Implementierung der rel_errorFunktion korrekt? Sollte ich nur überprüfen, ob abs(a-b)es weniger als winzig (= denormal) ist, und 0 zurückgeben? Oder sollte ich eine andere Kombination überprüfen, wie max(abs(a), abs(b))?

Ich möchte nur wissen, was der "richtige" Weg ist.

Ondřej Čertík
quelle

Antworten:

11

Sie können direkt für denormals überprüfen Verwendung isnormal()von math.h(C99 oder später POSIX.1 oder später). Wenn das Modul in Fortran ieee_arithmeticverfügbar ist, können Sie es verwenden ieee_is_normal(). Um die Fuzzy-Gleichheit genauer zu bestimmen, müssen Sie die Gleitkommadarstellung von Denormalen berücksichtigen und entscheiden, was Sie meinen, damit die Ergebnisse gut genug sind.

Um zu glauben, dass eines der beiden Ergebnisse korrekt ist, müssen Sie sicher sein, dass Sie in einem Zwischenschritt nicht zu viele Ziffern verloren haben. Das Rechnen mit Denormalen ist im Allgemeinen unzuverlässig und sollte vermieden werden, indem Ihr Algorithmus intern neu skaliert wird. Um sicherzustellen, dass Ihre interne Skalierung erfolgreich war, empfehle ich, Gleitkomma-Ausnahmen feenableexcept()mit C99 oder dem ieee_arithmeticModul in Fortran zu aktivieren .

Obwohl Ihre Anwendung das Signal abfangen kann, das bei Gleitkomma-Ausnahmen ausgelöst wird, setzen alle Kernel, die ich versucht habe, das Hardware-Flag zurück, sodass fetestexcept()kein nützliches Ergebnis zurückgegeben wird. Bei der Ausführung mit -fp_trapdrucken PETSc-Programme (standardmäßig) eine Stapelverfolgung, wenn ein Gleitkommafehler ausgelöst wird, identifizieren jedoch nicht die fehlerhafte Zeile. Wenn Sie in einem Debugger ausgeführt werden, behält der Debugger das Hardware-Flag bei und bricht den fehlerhaften Ausdruck ab. Sie können den genauen Grund überprüfen, indem Sie fetestexceptvom Debugger aus aufrufen, wobei das Ergebnis ein bitweises ODER der folgenden Flags ist (Werte können je nach Maschine variieren, siehe fenv.h; diese Werte gelten für x86-64 mit glibc).

  • FE_INVALID = 0x1
  • FE_DIVBYZERO = 0x4
  • FE_OVERFLOW = 0x8
  • FE_UNDERFLOW = 0x10
  • FE_INEXACT = 0x20
Jed Brown
quelle
Danke für die hervorragende Antwort. Der analytische Ausdruck, mit dem ich im asymptotischen Regime vergleiche, ist exp(log_gamma(m+0.5_dp) - (m+0.5_dp)*log(t)) / 2für m = 234, t = 2000. Es geht schnell auf Null, wenn ich zunehme m. Ich möchte nur sicherstellen, dass meine numerische Routine "korrekte" Zahlen (die Rückgabe von Null ist auch vollkommen in Ordnung) auf mindestens 12 signifikante Stellen zurückgibt. Wenn die Berechnung also eine Denormalzahl zurückgibt, ist sie einfach Null, und es sollte kein Problem geben. Daher muss nur die Vergleichsroutine robust sein.
Ondřej Čertík
5

Donald Knuth hat einen Vorschlag für einen Gleitkomma-Vergleichsalgorithmus in Band 2 "Seminumerische Algorithmen" von "The Art of Computer Programming". Es wurde in C von Th implementiert. Belding (siehe fcmp-Paket ) und ist in der GSL verfügbar .

GertVdE
quelle
2
Hier ist meine Fortran-Implementierung: gist.github.com/3776847 , beachten Sie, dass ich ohnehin explizit mit denormalen Zahlen umgehen muss. Ansonsten denke ich, dass es dem relativen Fehler ziemlich äquivalent ist. Der einzige Unterschied besteht darin, dass abs(a-b)/max(a, b) < epswir es tun abs(a-b)/2**exponent(max(a, b)) < eps, anstatt es zu tun , was die Mantisse max(a, b)so ziemlich einfach fallen lässt , also ist der Unterschied meiner Meinung nach vernachlässigbar.
Ondřej Čertík
5

Optimal gerundete denormalisierte Zahlen können tatsächlich einen hohen relativen Fehler aufweisen. (Das auf Null zu setzen, während es immer noch als relativer Fehler bezeichnet wird, ist irreführend.)

Nahe Null ist die Berechnung relativer Fehler jedoch bedeutungslos.

Daher sollten Sie, noch bevor Sie denormalisierte Zahlen erreichen, wahrscheinlich auf absolute Genauigkeit umschalten (nämlich die, die Sie in diesem Fall garantieren möchten).

yx|yx|absacc+relaccmax(|x|,|y|)

Dann wissen Benutzer Ihres Codes genau, wie genau sie wirklich sind.

Arnold Neumaier
quelle
Sind Sie sicher, dass es sinnlos ist, relative Fehler nahe Null zu berechnen? Ich denke, es ist nur dann bedeutungslos, wenn es zu einem Genauigkeitsverlust kommt (aus welchem ​​Grund auch immer). Wenn zum Beispiel aufgrund einiger numerischer Probleme (wie dem Subtrahieren von zwei großen Zahlen) ein Genauigkeitsverlust für x <1e-150 auftritt, haben Sie Recht. In meinem Fall scheinen die Zahlen jedoch bis auf Null genau zu sein, außer wenn sie die denormalen Zahlen treffen. Also in meinem Fall absacc = 1e-320 oder so und ich kann es einfach so überprüfen, abs(a-b) < tiny(1._dp)wie ich es oben mache.
Ondřej Čertík
@ OndřejČertík: In diesem Fall ersetzen Sie die 1e-150 durch 1e-300 oder was auch immer Sie überprüfen können. In jedem Fall machen Sie einen absoluten Fehler, der sehr nahe bei Null liegt, und Ihr Fehleranspruch sollte dies widerspiegeln, anstatt den relativen Fehler als Null zu deklarieren.
Arnold Neumaier
Aha. Ich kann überprüfen, ob alle für höhere Zahlen als funktionieren tiny(1._dp)=2.22507385850720138E-308(ich habe in meinem vorherigen Kommentar einen Fehler gemacht, es ist 2e-308, nicht 1e-320). Das ist also mein absoluter Fehler. Dann muss ich den relativen Fehler vergleichen. Ich verstehe deinen Standpunkt, ich denke du hast recht. Vielen Dank!
Ondřej Čertík
1
@ OndřejČertík: Um den zusätzlichen relativen Fehler bei absacc zu ermitteln, überwachen Sie das Maximum von . |yx|absaccmax(|x|,|y|)
Arnold Neumaier