Robuste Berechnung des Mittelwerts zweier Zahlen im Gleitkomma?

15

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.

becko
quelle
(+1) Interessante Frage, überraschend nicht trivial.
Kirill
1
In der Vergangenheit wurden Gleitkommawerte berechnet und für Zwischenergebnisse in einer präziseren Form gehalten. Wenn a + b (64-Bit-Doubles) ein 80-Bit-Zwischenergebnis ergibt und dies durch 2 geteilt wird, müssen Sie sich keine Gedanken über einen Überlauf machen. Ein Präzisionsverlust ist weniger offensichtlich.
JDługosz
Die Lösung hierfür scheint relativ einfach zu sein ( ich habe eine Antwort hinzugefügt ). Die Sache ist, ich bin ein Programmierer und kein Informatik-Experte, also was fehle ich, was diese Frage so viel schwieriger macht?
IQAndreas
Machen Sie sich keine Sorgen über die Kosten für Multiplikationen und Divisionen durch zwei. Ihr Compiler wird sie für Sie optimieren.
Federico Poloni

Antworten:

18

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, ob0xyx z y erfüllt. Der folgende Haskell-Codez=(x+y)/2xzy

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 .xfun(x,y)yx,y0xy

λ> 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)×2x

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 )?xy0xy

λ> 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 xx+(y/2-x/2)ySFloatSDouble

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.

Kirill
quelle
2
Warten Sie mal; Haben Sie behauptet, dass x + (y / 2-x / 2) eine gute Methode ist, wenn x <= y ist (unabhängig davon, ob x> = 0 ist oder nicht)? Scheint mir, dass das nicht richtig sein kann, da es im folgenden Fall die falsche Antwort gibt, wenn die Antwort genau darstellbar ist: x = -1, y = 1 + 2 ^ -52 (die kleinste darstellbare Zahl größer als 1), In diesem Fall lautet die Antwort 2 ^ -53. Bestätigung in Python: >>> x = -1.; y = 1.+2.**-52; print `2**-53`, `(x+y)/2.`, `x+(y/2.-x/2.)`
Don Hatch
2
x(x+y)/2yx,y(x+y)/2(x+y)/2
8

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.
Don Hatch
quelle
Mir ist klar, dass ich mit "denormalisiert" wirklich etwas anderes gemeint habe - dh Zahlen, die so nah beieinander liegen wie Zahlen, dh der Bereich von Zahlen, der ungefähr doppelt so groß ist wie der Bereich von denormalisierten Zahlen. dh die ersten 8 Ticks oder so im Diagramm unter en.wikipedia.org/wiki/Denormal_number . Der Punkt ist, die "ungeraden" von diesen sind die einzigen Zahlen, für die die Teilung durch zwei nicht genau ist. Ich muss diesen Teil der Antwort umformulieren, um dies zu verdeutlichen.
Don Hatch
fl(Öp(x,y))=Öp(x,y)(1+δ)|δ|ux/2+y/2(x+y)/2sind immer richtig gerundet, Über- / Unterlauf fehlt, es bleibt nur nichts Über- / Unterlauf zu zeigen, was einfach ist.
Kirill
@Kirill Ich bin ein bisschen verloren ... wo kommst du her? Ich denke auch nicht, dass es ganz richtig ist, dass "Divisionen durch 2 für nicht denormalen Zahlen genau sind" ... das ist das gleiche, worüber ich gestolpert bin, und es scheint ein bisschen umständlich zu sein, zu versuchen, es richtig zu machen. Die genaue Aussage ist eher wie "x / 2 ist genau, solange abs (x) mindestens doppelt so groß ist wie die größte subnormale Zahl" ... äh, peinlich!
Don Hatch
3

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)/2x/2+y/2binary64C[2-967,2970]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 )

njuffa
quelle
2

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.

IQAndreas
quelle
Ich versuche jetzt herauszufinden, wie dieselbe Antwort mit mehr als zwei Elementen funktioniert , wobei alle Variablen unter der größten Zahl bleiben und nur eine Teilungsoperation verwendet wird, um die Genauigkeit zu gewährleisten.
IQAndreas
@becko Yup, du würdest mindestens zweimal aufteilen. Auch das Beispiel, das Sie gaben, würde die Antwort falsch erscheinen lassen. Stellen Sie sich das Mittel von vor 2,4,9, es ist nicht dasselbe wie das Mittel von 3,9.
IQAndreas
Du hast recht, meine Rekursion war falsch. Ich bin mir nicht sicher, wie ich es im Moment beheben kann, ohne an Präzision zu verlieren.
Becko
Können Sie nachweisen, dass dies zu einem möglichst genauen Ergebnis führt? Das heißt, wenn xund ysind Gleitkomma, Ihre Berechnung erzeugt ein Gleitkomma am nächsten zu (x+y)/2?
Becko
1
Wird dies nicht ein Überlauf sein, wenn x, y die kleinste und größte auszudrückende Zahl sind?
Don Hatch
1

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.

leeroy
quelle
Dies ist der Brute-Force-Ansatz. Es funktioniert wahrscheinlich, aber ich habe nach einer Analyse gesucht, die keine höhere Präzision erfordert. Können Sie auch abschätzen, wie viel höhere Präzision zwischenzeitlich erforderlich ist? Lösche diese Antwort auf keinen Fall (+1), ich akzeptiere sie einfach nicht als Antwort.
Becko
1

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.

Roland Heide
quelle
0

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.
Keine Antwort
quelle