Solidne obliczenie średniej z dwóch liczb zmiennoprzecinkowych?

15

Pozwól x, ybyć dwóch liczb zmiennoprzecinkowych. Jaki jest właściwy sposób na obliczenie ich średniej?

Naiwny sposób (x+y)/2może doprowadzić do przepełnienia kiedy xi ysą zbyt duże. Myślę, że 0.5 * x + 0.5 * ymoże lepiej, ale wiąże się to z dwoma multiplikacjami (co może być nieefektywne) i nie jestem pewien, czy jest wystarczająco dobry. Czy jest lepszy sposób?

Innym pomysłem, z którym bawiłem się, jest to, (y/2)(1 + x/y)czy x<=y. Ale znowu nie jestem pewien, jak to przeanalizować i udowodnić, że spełnia moje wymagania.

Ponadto potrzebuję gwarancji, że obliczoną średnią będzie >= min(x,y)i <= max(x,y). Jak wskazano w odpowiedzi Dona Hatcha , może lepszym sposobem postawienia tego pytania jest: Jaka jest implementacja średniej z dwóch liczb, która zawsze daje najbardziej możliwy dokładny wynik? To znaczy, jeśli xi yczy liczby zmiennoprzecinkowe są najbliżej, jak obliczyć liczbę zmiennoprzecinkową (x+y)/2? W takim przypadku obliczona średnia jest automatycznie >= min(x,y)i <= max(x,y). Szczegółowe informacje można znaleźć w odpowiedzi Dona Hatcha .

Uwaga: Moim priorytetem jest solidna dokładność. Wydajność jest zbywalna. Jeśli jednak istnieje wiele niezawodnych i dokładnych algorytmów, wybrałbym najbardziej wydajny.

becko
źródło
(+1) Ciekawe pytanie, zaskakująco nietrywialne.
Kirill,
1
W przeszłości wartości zmiennoprzecinkowe były obliczane i utrzymywane w formie o wyższej precyzji dla wyników pośrednich. Jeśli a + b (64-bitowe podwojenie) daje 80-bitowy wynik pośredni i to jest podzielone przez 2, nie musisz się martwić o przepełnienie. Utrata precyzji jest mniej oczywista.
JDługosz
Rozwiązanie tego wydaje się stosunkowo proste ( dodałem odpowiedź ). Chodzi o to, że jestem programistą, a nie ekspertem w dziedzinie informatyki, więc czego mi brakuje, co sprawia, że ​​to pytanie jest o wiele trudniejsze?
IQAndreas
Nie martw się kosztem mnożenia i dzielenia przez dwa; Twój kompilator zoptymalizuje je dla Ciebie.
Federico Poloni

Odpowiedzi:

18

Myślę, że dokładność i stabilność algorytmów numerycznych Highama dotyczy tego, w jaki sposób można analizować tego rodzaju problemy. Patrz rozdział 2, zwłaszcza ćwiczenie 2.8.

W tej odpowiedzi chciałbym wskazać coś, co tak naprawdę nie zostało poruszone w książce Highama (wydaje się, że nie jest to zbyt powszechnie znane). Jeśli jesteś zainteresowany udowodnieniem właściwości prostych algorytmów numerycznych takich jak te, możesz skorzystać z mocy nowoczesnych solverów SMT ( Teorie satysfakcji modulo ), takich jak z3 , używając pakietu takiego jak sbv w Haskell. Jest to nieco łatwiejsze niż używanie ołówka i papieru.

Załóżmy, że podano mi i chciałbym wiedzieć, czy spełnia . Poniższy kod Haskellz = ( x + y ) / 2 x z y0xyz=(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

pozwoli mi to zrobić automatycznie . Tutaj test1 funjest założenie , że dla wszystkich skończonej pływaki z .x , y 0 x yxfun(x,y)yx,y0xy

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

To się przelewa. Załóżmy, że teraz biorę inną formułę: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

Nie działa (ze względu na stopniowe niedopełnienie: , co może być nieintuicyjne, ponieważ cała arytmetyka ma wartość base-2).(x/2)×2x

Teraz spróbuj :z=x+(yx)/2

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

Pracuje! Q.E.D.To dowód , że test1własność zachodzi dla wszystkich pływaków, jak zdefiniowano powyżej.

Co z tym samym, ale ograniczonym do (zamiast )?xy0xy

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

Okej, więc jeśli przepełni, co powiesz na ?yxz=x+(y/2x/2)

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

Wygląda więc na to, że spośród wzorów, które tu wypróbowałem, wydaje się działać (również z dowodem). Metoda solvera SMT wydaje mi się o wiele szybszym sposobem odpowiedzi na podejrzenia dotyczące prostych wzorów zmiennoprzecinkowych niż analizowanie błędów zmiennoprzecinkowych ołówkiem i papierem.x+(y/2x/2)

Wreszcie cel dokładności i stabilności często stoi w sprzeczności z celem wydajności. Jeśli chodzi o wydajność, tak naprawdę nie widzę, jak możesz sobie radzić lepiej niż , zwłaszcza, że ​​kompilator nadal będzie cię ciężko tłumaczyć, tłumacząc to na instrukcje maszynowe.(x+y)/2

PS Wszystko to z arytmetyką zmiennoprzecinkową IEEE754 o pojedynczej precyzji. Sprawdziłem z podwójnej precyzji arytmetyki (wymienić z ) i działa zbyt.xx+(y/2x/2)ySFloatSDouble

PPS Jedną z rzeczy, o których należy pamiętać przy implementacji tego w kodzie jest to, że flagi kompilatora takie jak -ffast-math(niektóre formy takich flag są czasami domyślnie włączone w niektórych popularnych kompilatorach) nie spowodują arytmetyki IEEE754, co unieważni powyższe dowody. Jeśli używasz flag, które umożliwiają np. Optymalizacje dodawania skojarzonego, nie ma sensu robić niczego innego niż .(x+y)/2

PPPS Dałam się trochę ponieść spojrzeniu tylko na proste wyrażenia algebraiczne bez warunków warunkowych. Don Hatch „s formuła jest ściśle lepiej.

Cyryl
źródło
2
Czekaj; czy twierdziłeś, że jeśli x <= y (niezależnie od tego, czy x> = 0 czy nie), to x + (y / 2-x / 2) jest dobrym sposobem na zrobienie tego? Wydaje mi się, że to nieprawda, ponieważ daje złą odpowiedź w następującym przypadku, gdy odpowiedź jest dokładnie reprezentowalna: x = -1, y = 1 + 2 ^ -52 (najmniejsza reprezentowalna liczba większa niż 1), w takim przypadku odpowiedź to 2 ^ -53. Potwierdzenie w 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

Po pierwsze, zauważ, że jeśli masz metodę, która daje najdokładniejszą odpowiedź we wszystkich przypadkach, spełni ona wymagany warunek. (Należy pamiętać, że mówię najdokładniejszą odpowiedź zamiast z najdokładniejszej odpowiedzi, ponieważ nie może być dwóch zwycięzców.) Dowód: Jeśli, przeciwnie, trzeba dokładnej-as-możliwa odpowiedź, która ma nie spełniają wymaganego warunku, że oznacza albo (w którym przypadku jest lepsza odpowiedź, sprzeczność), albo (w którym przypadku jest lepsza odpowiedź, sprzeczność).answer<min(x,y)<=max(x,y)min(x,y)min(x,y)<=max(x,y)<answermax(x,y)

Myślę więc, że to oznacza, że ​​twoje pytanie sprowadza się do znalezienia najdokładniejszej możliwej odpowiedzi. Zakładając, że arytmetyka IEEE754 jest w toku, proponuję:

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

Mój argument, że daje to najdokładniejszą odpowiedź, jest dość żmudną analizą przypadku. Tutaj idzie:

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

    • Podklucz ani x, ani y nie są denormalizowane: w tym przypadku obliczona odpowiedź x/2.+y/2.manipuluje tymi samymi mantysami, a zatem daje dokładnie taką samą odpowiedź, jak obliczenie (x+y)/2dałoby, gdybyśmy przyjęli rozszerzone wykładniki, aby zapobiec przepełnieniu. Ta odpowiedź może zależeć od trybu zaokrąglania, ale w każdym przypadku IEEE754 gwarantuje, że jest to najlepsza możliwa odpowiedź (z faktu, że obliczone x+ygwarantuje najlepsze przybliżenie do matematyki x + y, a podział przez 2 jest w tym dokładny walizka).
    • Podtekst x jest zdenormalizowany (i tak 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.

    • Podsekcja y jest zdenormalizowana (i tak abs(x)>=1): analogicznie.

  • Sprawa max(abs(x),abs(y)) < 1.:
    • Podliczona obliczona x+yjest albo niez denormalizowana, albo zdenormalizowana, a „- nawet”: Chociaż obliczona x+ymoże nie być dokładna, IEEE754 gwarantuje, że jest najlepszym możliwym przybliżeniem do matematyki x + y. W tym przypadku kolejny podział przez 2 w wyrażeniu (x+y)/2.jest dokładny, więc obliczona odpowiedź (x+y)/2.jest najlepszym możliwym przybliżeniem do matematyki (x + y) / 2.
    • Podklucz obliczonego x+yjest zdenormalizowany i „nieparzysty”: W tym przypadku dokładnie jeden z x, y również musi być zdenormalizowany - i - „nieparzysty”, co oznacza, że ​​drugi z x, y jest zdenormalizowany znakiem przeciwnym, a zatem obliczone x+yjest dokładnie matematyczne x + y, a zatem obliczone (x+y)/2.jest gwarantowane przez IEEE754, aby być najlepszym możliwym przybliżeniem do matematyki (x + y) / 2.
Don Hatch
źródło
Zdaję sobie sprawę, kiedy powiedziałem „zdenormalizowany”, naprawdę miałem na myśli coś innego - to znaczy liczby, które są tak blisko siebie, jak to tylko możliwe, tj. Zakres liczb, który jest około dwa razy większy niż zakres liczb zdenormalizowanych, tj. pierwszych 8 tyknięć na diagramie w en.wikipedia.org/wiki/Denormal_number . Chodzi o to, że „nieparzyste” z nich są jedynymi liczbami, dla których dzielenie ich przez dwa nie jest dokładne. Muszę ponownie sformułować tę część odpowiedzi, aby było to jasne.
Don Hatch,
fl(op(x,y))=op(x,y)(1+δ)|δ|ux/2+y/2(x+y)/2są zawsze poprawnie zaokrąglone, brak przepełnienia / niedopełnienia, wszystko, co pozostaje, to nie pokazywanie niczego, co jest przepełnieniem / niedopełnieniem, co jest łatwe.
Kirill,
@Kirill Jestem trochę zagubiony ... skąd pochodzisz? Nie sądzę też, żeby to prawda, że ​​„podział na 2 jest dokładny dla liczb nienormalnych” ... to jest to samo, o co się potknąłem i wydaje się to trochę dziwne, aby spróbować to naprawić. Dokładne stwierdzenie jest bardziej podobne do tego, że „x / 2 jest dokładne, o ile abs (x) jest co najmniej dwa razy większa od największej liczby nienormalnej” ... argh, niezręczne!
Don Hatch,
3

W przypadku binarnych formatów zmiennoprzecinkowych IEEE-754, których przykładem jest binary64obliczenie (podwójnej precyzji), S. Boldo formalnie udowodnił, że przedstawiony poniżej prosty algorytm zapewnia poprawnie zaokrągloną średnią.

Sylvie Boldo, „Formalna weryfikacja programów obliczających średnią zmiennoprzecinkową”. Na międzynarodowej konferencji na temat formalnych metod inżynieryjnych , s. 17–32. Springer, Cham, 2015. ( projekt online )

(x+y)/2x/2+y/2binary64C[2967,2970]C tak aby zapewnić najlepszą wydajność dla konkretnego przypadku użycia.

Daje to następujący przykładowy ISO-C99kod:

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);
}

W ostatnich pracach uzupełniających S. Boldo i współautorzy pokazali, jak osiągnąć najlepsze możliwe wyniki dla dziesiętnych formatów zmiennoprzecinkowych IEEE-754, wykorzystując operacje fuzji wielokrotnego dodawania (FMA) i dobrze znaną precyzję podwajanie bloku konstrukcyjnego (TwoSum):

Sylvie Boldo, Florian Faissole i Vincent Tourneur, „Formalnie ustalony algorytm obliczania poprawnej średniej liczb dziesiętnych zmiennoprzecinkowych”. W 25. sympozjum IEEE na temat arytmetyki komputerowej (ARITH 25) , czerwiec 2018 r., S. 69–75. ( projekt online )

njuffa
źródło
2

Chociaż nie może być super-wydajny wydajność mądry, istnieje bardzo prosty sposób (1) upewnij się, że żaden z tych liczb jest większa niż którakolwiek xlub y(bez przepełnienia) i (2) utrzymać zmiennoprzecinkowych jako „dokładne”, jak możliwe (i (3) , jako dodatkowy bonus, nawet jeśli stosowane jest odejmowanie, żadne wartości nie będą nigdy przechowywane jako liczby ujemne.

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

W rzeczywistości, jeśli naprawdę chcesz uzyskać dokładność, nie musisz nawet dokonywać podziału na miejscu; po prostu zwróć wartości min(x, y)i differencektórych możesz użyć, aby uprościć logicznie lub manipulować później.

IQAndreas
źródło
Teraz próbuję wymyślić, jak sprawić, by ta sama odpowiedź działała z więcej niż dwoma elementami , jednocześnie utrzymując wszystkie zmienne na poziomie niższym od największej z liczb i używając tylko jednej operacji dzielenia, aby zachować dokładność.
IQAndreas
@becko Tak, dzieliłbyś się co najmniej dwa razy. Podany przykład sprawiłby, że odpowiedź byłaby błędna. Wyobraź sobie średnią 2,4,9, to nie to samo co średnia 3,9.
IQAndreas
Masz rację, moja rekurencja była błędna. Nie jestem pewien, jak to naprawić teraz, bez utraty precyzji.
becko
Czy możesz udowodnić, że daje to jak najdokładniejszy wynik? To znaczy, jeśli xi yczy zmiennoprzecinkowe, twoje obliczenia dają zmiennoprzecinkowy najbliższy (x+y)/2?
becko
1
Czy to nie przepełni się, gdy x, y są najmniejszymi i największymi możliwymi do wyrażenia liczbami?
Don Hatch,
1

Konwertuj na wyższą precyzję, dodaj tam wartości i przekonwertuj z powrotem.

Wyższa precyzja nie powinna powodować przepełnienia, a jeśli oba są w prawidłowym zakresie zmiennoprzecinkowym, obliczona liczba również powinna być w środku.

I powinna znajdować się pomiędzy nimi, w najgorszym przypadku tylko połowa większej liczby, jeśli prewencja nie jest wystarczająca.

leeroy
źródło
To podejście brutalnej siły. Prawdopodobnie działa, ale szukałem analizy, która nie wymagałaby pośredniej wyższej precyzji. Czy potrafisz także oszacować, jaka jest wymagana pośrednia wyższa precyzja? W każdym razie nie usuwaj tej odpowiedzi (+1), po prostu nie zaakceptuję jej jako odpowiedzi.
becko
1

Teoretycznie x/2można go obliczyć, odejmując 1 od mantysy.

Jednak faktyczne implementowanie takich operacji bitowych niekoniecznie jest proste, szczególnie jeśli nie znasz formatu liczb zmiennoprzecinkowych.

Jeśli możesz to zrobić, cała operacja zostaje zredukowana do 3 dodawania / odejmowania, co powinno być znaczącą poprawą.

Roland Heath
źródło
0

Myślałem w tym samym kierunku, co @Roland Heath, ale nie mogę jeszcze komentować, oto moje zdanie:

x/2można obliczyć, odejmując 1 od wykładnika (nie mantysy, odejmując 1 od mantysy odejmując 2^(value_of_exponent-length_of_mantissa)od wartości ogólnej).

Załóżmy, że bez ograniczenia ogólnego przypadku x < y. (Jeśli x > y, ponownie oznacz zmienne. Jeśli x = y, (x+y) / 2jest banalne.)

  • Przekształć (x+y) / 2w x/2 + y/2, który można wykonać przez dwie odejmowanie liczb całkowitych (jedną od wykładnika)
    • Jednak wykładnik ma dolną granicę w zależności od reprezentacji. Jeśli wykładnik jest już minimalny przed odjęciem 1, ta metoda będzie wymagać specjalnej obsługi przypadków. Minimalny wykładnik włączony xsprawi, że będzie on x/2mniejszy niż reprezentowalny (zakładając, że mantysa jest reprezentowana z ukrytym prowadzeniem 1).
    • Zamiast odejmować 1 od wykładnika x, xprzesuń mantysę w prawo o jeden (i dodaj ukrytą wiodącą 1, jeśli istnieje).
    • Odejmij 1 od wykładnika y, jeśli nie jest minimalny. Jeśli jest minimalna (y jest większe niż x, ze względu na mantysę), przesuń mantysę w prawo o jeden (dodaj ukryte prowadzenie 1, jeśli istnieje).
    • Przesuń nową mantysę xw prawo zgodnie z wykładnikiem y.
    • Wykonaj dodawanie liczb całkowitych na mantissae, chyba że mantysa xzostała całkowicie przesunięta. Jeśli oba wykładniki były minimalne, wiodące przepełniłyby się, co jest w porządku, ponieważ to przepełnienie powinno stać się ponownie domyślnym wiodącym.
  • i jeden dodatek zmiennoprzecinkowy.
    • Nie mogę tutaj wymyślić żadnego specjalnego przypadku; z wyjątkiem zaokrąglania, które ma również zastosowanie do przesunięć opisanych powyżej.
Brak odpowiedzi
źródło