Robuste Berechnung des Mittelwerts zweier Zahlen im Gleitkomma?

Lassen Sie x, ywerden zwei Gleitkommazahlen. Wie kann man den Mittelwert richtig berechnen?

Der naive Weg (x+y)/2kann zu Überläufen führen, wenn xund ysind zu groß. Ich denke 0.5 * x + 0.5 * yvielleicht besser, aber es geht um zwei Multiplikationen (was vielleicht ineffizient ist), und ich bin nicht sicher, ob es gut genug ist. Gibt es einen besseren Weg?

Eine andere Idee, mit der ich gespielt habe, ist (y/2)(1 + x/y)wenn x<=y. Aber auch hier bin ich mir nicht sicher, wie ich das analysieren und beweisen soll, dass es meinen Anforderungen entspricht.

Außerdem brauche ich eine Garantie, dass der berechnete Mittelwert >= min(x,y)und ist <= max(x,y). Wie in der Antwort von Don Hatch ausgeführt , lautet eine bessere Möglichkeit, diese Frage zu stellen: Was ist eine Implementierung des Mittelwerts aus zwei Zahlen, die immer das genaueste Ergebnis liefert? Das heißt, wenn xund ysind Gleitkommazahlen, wie wird die Gleitkommazahl berechnet, die am nächsten ist (x+y)/2? In diesem Fall ist der berechnete Mittelwert automatisch >= min(x,y)und <= max(x,y). Siehe Don Hatchs Antwort für Details.

Hinweis: Meine Priorität ist robuste Genauigkeit. Effizienz ist entbehrlich. Wenn es jedoch viele robuste und genaue Algorithmen gibt, würde ich die effizienteste auswählen.

Ich denke, Highams Genauigkeit und Stabilität numerischer Algorithmen befassen sich damit, wie man diese Art von Problemen analysieren kann. Siehe Kapitel 2, insbesondere Übung 2.8.

In dieser Antwort möchte ich auf etwas hinweisen, das in Highams Buch nicht wirklich angesprochen wird (es scheint im Übrigen nicht sehr bekannt zu sein). Wenn Sie interessiert sind zu beweisen Eigenschaften von einfachen numerischen Algorithmen wie diese, können Sie die Macht der modernen SMT - Solver (verwenden Satisfiability Modulo Theories ), wie z3 , ein Paket wie die Verwendung von SBV in Haskell. Das ist etwas einfacher als mit Bleistift und Papier.

Angenommen, mir wird , und ich würde gerne wissen, ob$0\leq x\leq y$x ≤ z ≤ y erfüllt. Der folgende Haskell-Code$z=(x+y)/2$$x\leq z\leq y$

import Data.SBV

test1 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test1 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ 0 .<= x &&& x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

test2 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test2 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

werde mich das automatisch machen lassen . Hier test1 funist der Satz, dass für alle finiten floats x , y mit 0 ≤ x ≤ y .$x \leq \mathit{fun}(x,y) \leq y$$x,y$$0\leq x\leq y$

λ> prove $ test1 (\x y -> (x + y) / 2)
Falsifiable. Counter-example:
  x = 2.3089316e36 :: Float
  y = 3.379786e38 :: Float

Es läuft über. Angenommen, ich nehme jetzt Ihre andere Formel: $z=x/2+y/2$

λ> prove $ test1 (\x y -> x/2 + y/2)
Falsifiable. Counter-example:
  x = 2.3509886e-38 :: Float
  y = 2.3509886e-38 :: Float

Funktioniert nicht (aufgrund eines allmählichen Unterlaufs: , was möglicherweise nicht intuitiv ist, da alle Arithmetik zur Basis 2 gehört).$(x/2)\times2 \neq x$

Versuchen Sie nun :$z=x + (y-x)/2$

λ> prove $ test1 (\x y -> x + (y-x)/2)
Q.E.D.

Funktioniert! Dies Q.E.D.ist ein Beweis dafür, dass die test1Eigenschaft für alle oben definierten Floats gilt.

Was ist mit dem gleichen, aber auf (anstelle von 0 ≤ x ≤ y )?$x\leq y$$0\leq x\leq y$

λ> prove $ test2 (\x y -> x + (y-x)/2)
Falsifiable. Counter-example:
  x = -3.1300826e34 :: Float
  y = 3.402721e38 :: Float

Okay, wenn überläuft, wie wäre es dann mit z $y-x$ ?$z = x + (y/2-x/2)$

λ> prove $ test2 (\x y -> x + (y/2 - x/2))
Q.E.D.

So scheint es, dass unter den Formeln, die ich hier ausprobiert habe, zu funktionieren scheint (auch mit einem Beweis). Der SMT-Solver-Ansatz scheint mir eine viel schnellere Möglichkeit zu sein, Verdacht auf einfache Fließkommaformeln zu äußern, als die Fließkomma-Fehleranalyse mit Bleistift und Papier durchzuführen.$x + (y/2 - x/2)$

Schließlich steht das Ziel der Genauigkeit und Stabilität häufig im Widerspruch zum Ziel der Leistung. Was die Leistung angeht, sehe ich nicht wirklich, wie Sie es besser machen können als , zumal der Compiler immer noch die Mühe macht, dies in Maschinenanweisungen für Sie zu übersetzen.$(x+y)/2$

PS Dies ist alles mit IEEE754-Gleitkomma-Arithmetik mit einfacher Genauigkeit. Ich habe $x \leq x + (y/2-x/2) \leq y$SFloatSDouble

PPS-ffast-math$(x+y)/2$

PPPS Ich wurde ein wenig mitgerissen, als ich mir nur einfache algebraische Ausdrücke ohne Bedingungen ansah . Die Formel von Don Hatch ist streng besser.

Beachten Sie zunächst, dass eine Methode, die in allen Fällen die genaueste Antwort liefert, Ihre erforderliche Bedingung erfüllt. (Beachten Sie, dass ich sage , eine genaueste Antwort eher als die genaueste Antwort, da es kann zwei Gewinner sein.) Beweis: Wenn, im Gegenteil, Sie haben eine genaue as mögliche Antwort , die nicht nicht die erforderliche Bedingung erfüllen, dass bedeutet entweder answer<min(x,y)<=max(x,y)(in welchem ​​Fall min(x,y)ist eine bessere Antwort, ein Widerspruch) oder min(x,y)<=max(x,y)<answer(in welchem ​​Fall max(x,y)ist eine bessere Antwort, ein Widerspruch).

Ich denke, das bedeutet, dass Ihre Frage darauf hinausläuft, eine möglichst genaue Antwort zu finden. Unter der Annahme von IEEE754-Arithmetik schlage ich Folgendes vor:

if max(abs(x),abs(y)) >= 1.:
    return x/2. + y/2.
else:
    return (x+y)/2.

Mein Argument, dass dies die genaueste Antwort liefert, ist eine etwas langwierige Fallanalyse. Hier geht:

• Fall max(abs(x),abs(y)) >= 1.:

  • Unterfall weder x noch y ist denormalisiert: In diesem Fall x/2.+y/2.manipuliert die berechnete Antwort die gleichen Mantissen und gibt daher genau die gleiche Antwort wie die Berechnung von(x+y)/2 würde ergeben, wenn wir erweiterte Exponenten würden, um einen Überlauf zu verhindern. Diese Antwort kann vom Rundungsmodus abhängen, aber in jedem Fall garantiert IEEE754, dass es sich um eine bestmögliche Antwort handelt (aus der Tatsache, dass die berechnete x+yNäherung an das mathematische x + y garantiert ist und die Division durch 2 genau ist) Fall).

  • Der Unterfall x ist denormalisiert (und so abs(y)>=1):

    answer = x/2. + y/2. = y/2. since abs(x/2.) is so tiny compared to abs(y/2.) = the exact mathematical value of y/2 = a best possible answer.

  • Der Unterfall y ist denormalisiert (und so abs(x)>=1): analog.

• Fall max(abs(x),abs(y)) < 1. :
  • Der berechnete Teilfall x+yist entweder nicht-denormalisiert oder denormalisiert-und- "gerade": Obwohl der berechnete x+yWert möglicherweise nicht genau ist, wird durch IEEE754 eine bestmögliche Annäherung an das mathematische x + y garantiert. In diesem Fall ist die nachfolgende Division durch 2 im Ausdruck (x+y)/2.genau, sodass die berechnete Antwort (x+y)/2.eine bestmögliche Annäherung an das mathematische (x + y) / 2 darstellt.
  • Unterfall des berechneten x+ywird denormiert und „ungerade“: In diesem Fall genau ein von x, y muss auch denormalisierte-und- „odd“ sein , welches die andere von X bedeutet, y mit dem entgegengesetzten Vorzeichen denormalisiert ist, und so die berechnete x+yIST genau das mathematische x + y, und so (x+y)/2.wird durch IEEE754 garantiert , dass das berechnete eine bestmögliche Annäherung an das mathematische (x + y) / 2 ist.

Für binäre Gleitkommaformate nach IEEE-754, am Beispiel von binary64 (doppelte Genauigkeits-) Berechnung , hat S. Boldo formal bewiesen, dass der unten gezeigte einfache Algorithmus den korrekt gerundeten Durchschnitt liefert.

Sylvie Boldo, "Formale Überprüfung von Programmen, die den Gleitkomma-Durchschnitt berechnen." In International Conference on Formal Engineering Methods , S. 17-32. Springer, Cham, 2015. ( Entwurf online )

$(x+y)/2$$x/2 + y/2$binary64$C \in [2^{-967}, 2^{970}]$$C$ um die beste Leistung für einen bestimmten Anwendungsfall bereitzustellen.

Dies ergibt den folgenden beispielhaften ISO-C99Code:

double average (double x, double y) 
{
    const double C = 1; /* 0x1p-967 <= C <= 0x1p970 */
    return (C <= fabs (x)) ? (x / 2 + y / 2) : ((x + y) / 2);
}

In den jüngsten Nacharbeiten haben S. Boldo und Mitautoren gezeigt, wie die bestmöglichen Ergebnisse für die IEEE-754-Dezimal-Gleitkommaformate erzielt werden können, indem FMA-Operationen (Fused Multiply Add) und eine bekannte Präzisionsmethode verwendet werden. Baustein verdoppeln (TwoSum):

Sylvie Boldo, Florian Faissole und Vincent Tourneur, "Ein formell erprobter Algorithmus zur Berechnung des korrekten Durchschnitts von Gleitkommazahlen." Im 25. IEEE-Symposium für Computerarithmetik (ARITH 25) , Juni 2018, S. 69-75. ( Entwurf online )

Obwohl dies in Bezug auf die Leistung möglicherweise nicht besonders effizient ist, gibt es eine sehr einfache Möglichkeit, um (1) sicherzustellen, dass keine der Zahlen größer als entweder xoder y(keine Überläufe) ist, und (2) den Gleitkommawert so genau wie möglich zu halten möglich (und (3) , als zusätzlicher Bonus, obwohl Subtraktion verwendet wird, werden niemals Werte als negative Zahlen gespeichert.

float difference = max(x, y) - min(x, y);
return min(x, y) + (difference / 2.0);

In der Tat, wenn Sie wirklich Genauigkeit anstreben möchten, müssen Sie die Teilung nicht einmal an Ort und Stelle durchführen. Geben Sie einfach die Werte von min(x, y)und zurück, differencedie Sie zur Vereinfachung oder späteren Bearbeitung verwenden können.

In höhere Genauigkeit konvertieren, dort die Werte addieren und zurückkonvertieren.

Bei der höheren Genauigkeit sollte es keinen Überlauf geben, und wenn sich beide im gültigen Gleitkommabereich befinden, sollte die berechnete Zahl auch innerhalb liegen.

Und es sollte dazwischen liegen, im schlimmsten Fall nur die Hälfte der größeren Zahl, wenn die Präzision nicht ausreicht.

Theoretisch, x/2 kann durch Subtrahieren von 1 von der Mantisse berechnet werden.

Die tatsächliche Implementierung solcher bitweisen Operationen ist jedoch nicht unbedingt einfach, insbesondere wenn Sie das Format Ihrer Gleitkommazahlen nicht kennen.

Wenn Sie dies tun können, wird die gesamte Operation auf 3 Additionen / Subtraktionen reduziert, was eine signifikante Verbesserung darstellen sollte.

Ich habe nach dem Vorbild von @Roland Heath gedacht, kann mich aber noch nicht dazu äußern.

x/2kann durch Subtrahieren von 1 vom Exponenten berechnet werden (nicht von der Mantisse, Subtrahieren von 1 von der Mantisse ist Subtrahieren2^(value_of_exponent-length_of_mantissa) vom Gesamtwert ).

Nehmen wir ohne Einschränkung des allgemeinen Falls an x < y. (Wenn x > y, benennen Sie die Variablen neu. Wenn x = y, (x+y) / 2ist das trivial.)

• Verwandle dich (x+y) / 2inx/2 + y/2 , was durch zwei ganzzahlige Subtraktionen (durch eine vom Exponenten) durchgeführt werden kann
  • Abhängig von Ihrer Darstellung gibt es jedoch eine Untergrenze für den Exponenten. Wenn Ihr Exponent vor dem Subtrahieren von 1 bereits minimal ist, erfordert diese Methode eine spezielle Fallbehandlung. Ein minimaler Exponent xwird auf machenx/2 kleiner als darstellbar (vorausgesetzt, die Mantisse wird mit einer impliziten führenden 1 dargestellt).
  • Anstatt 1 vom Exponenten von zu subtrahieren x, verschieben Siex die Mantisse um eins nach rechts (und addieren Sie gegebenenfalls die implizite führende 1).
  • Subtrahiere 1 vom Exponenten von y, wenn es nicht minimal ist. Wenn es minimal ist (y ist wegen der Mantisse größer als x), verschieben Sie die Mantisse um eins nach rechts (fügen Sie gegebenenfalls implizite führende 1 hinzu).
  • Verschieben Sie die neue Mantisse von xnach rechts entsprechend dem Exponenten von y.
  • Führen Sie eine Ganzzahladdition an der Mantisse durch, es sei denn, die Mantisse von xwurde vollständig verschoben. Wenn beide Exponenten minimal wären, würden die führenden überlaufen, was in Ordnung ist, da dieser Überlauf wieder zu einem impliziten führenden werden soll.
• und eine Gleitkommazugabe.
  • Ich kann mir hier keinen besonderen Fall vorstellen. mit Ausnahme der Rundung, die auch für die oben beschriebene Verschiebung gilt.