Cálculo robusto de la media de dos números en coma flotante?

15

Dejado x, yser dos números de punto flotante. ¿Cuál es la forma correcta de calcular su media?

La forma ingenua (x+y)/2puede dar lugar a desbordamientos cuando xy yson demasiado grandes. Creo que 0.5 * x + 0.5 * ytal vez sea mejor, pero implica dos multiplicaciones (que tal vez sea ineficiente), y no estoy seguro de si es lo suficientemente bueno. ¿Hay una mejor manera?

Otra idea con la que he estado jugando es (y/2)(1 + x/y)si x<=y. Pero, una vez más, no estoy seguro de cómo analizar esto y demostrar que cumple con mis requisitos.

Además, necesito una garantía de que la media calculada será >= min(x,y)y <= max(x,y). Como se señaló en la respuesta de Don Hatch , quizás una mejor manera de plantear esta pregunta es: ¿Cuál es una implementación de la media de dos números que siempre da el resultado más exacto posible? Es decir, si xy yson números de coma flotante, ¿cómo calcular el número de coma flotante más cercano (x+y)/2? En este caso, la media calculada es automáticamente >= min(x,y)y <= max(x,y). Vea la respuesta de Don Hatch para más detalles.

Nota: Mi prioridad es la precisión robusta. La eficiencia es prescindible. Sin embargo, si hay muchos algoritmos robustos y precisos, elegiría el más eficiente.

llamar
fuente
(+1) Pregunta interesante, sorprendentemente no trivial.
Kirill
1
En el pasado, los valores de coma flotante se calculaban y mantenían en una forma de mayor precisión para obtener resultados intermedios. Si a + b (dobles de 64 bits) produce un resultado intermedio de 80 bits y esto es lo que se divide por 2, no tiene que preocuparse por el desbordamiento. La pérdida de precisión es menos obvia.
JDługosz
La solución a esto parece relativamente simple ( agregué una respuesta ). La cuestión es que soy un programador y no un experto en ciencias de la computación, entonces, ¿qué me estoy perdiendo que hace que esta pregunta sea mucho más difícil?
IQAndreas
No se preocupe por el costo de multiplicaciones y divisiones por dos; su compilador los optimizará por usted.
Federico Poloni

Respuestas:

18

Creo que la precisión y la estabilidad de los algoritmos numéricos de Higham aborda cómo uno puede analizar este tipo de problemas. Vea el Capítulo 2, especialmente el ejercicio 2.8.

En esta respuesta, me gustaría señalar algo que realmente no se aborda en el libro de Higham (no parece ser muy conocido, de hecho). Si está interesado en probar las propiedades de algoritmos numéricos simples como estos, puede usar el poder de los solucionadores SMT modernos ( Teorías del módulo de satisfacción ), como z3 , usando un paquete como sbv en Haskell. Esto es algo más fácil que usar lápiz y papel.

Supongamos que se me da que , y me gustaría saber si z = ( x + y ) / 2 satisface x z y . El siguiente código de Haskell0 0Xyz=(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

me dejará hacer esto automáticamente . Aquí test1 funestá la proposición de que para todos los flotadores finitos x , y con 0 x y .XFtunorte(X,y)yX,y0 0Xy

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

Se desborda. Supongamos que ahora tomo su otra fórmula: 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

No funciona (debido al flujo inferior gradual: , lo que podría no ser intuitivo debido a que toda la aritmética es base-2).(X/ /2)×2X

Ahora intente :z=X+(y-X)/ /2

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

¡Trabajos! El Q.E.D.es una prueba de que la test1propiedad es válida para todas las carrozas como se definió anteriormente.

¿Qué pasa con lo mismo, pero restringido a (en lugar de 0 x y )?Xy0 0Xy

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

Bien, entonces si desborda, ¿qué tal z = x + ( y / 2 - x / 2 ) ?y-Xz=X+(y/ /2-X/ /2)

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

Entonces parece que entre las fórmulas que he probado aquí, parece funcionar (con una prueba, también). El enfoque del solucionador SMT me parece una forma mucho más rápida de responder a las sospechas sobre fórmulas simples de punto flotante que pasar por un análisis de error de punto flotante con lápiz y papel.X+(y/ /2-X/ /2)

Finalmente, el objetivo de precisión y estabilidad a menudo está en desacuerdo con el objetivo de rendimiento. Para el rendimiento, realmente no veo cómo puede hacerlo mejor que , especialmente porque el compilador aún hará el trabajo pesado de traducir esto en instrucciones de máquina para usted.(X+y)/ /2

XX+(y/ /2-X/ /2)ySFloatSDouble

-ffast-math(X+y)/ /2

PPPS Me dejé llevar un poco mirando solo expresiones algebraicas simples sin condicionales. La fórmula de Don Hatch es estrictamente mejor.

Kirill
fuente
2
Espere; ¿Afirmó que si x <= y (independientemente de si x> = 0 o no), entonces x + (y / 2-x / 2) es una buena manera de hacerlo? Me parece que no puede ser correcto, ya que da la respuesta incorrecta en el siguiente caso cuando la respuesta es exactamente representable: x = -1, y = 1 + 2 ^ -52 (el número representable más pequeño mayor que 1), en cuyo caso la respuesta es 2 ^ -53. Confirmación en 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

Primero, observe que si tiene un método que da una respuesta más precisa en todos los casos, entonces satisfará su condición requerida. (Tenga en cuenta que digo una respuesta más precisa en lugar de la respuesta más precisa, ya que puede haber dos ganadores). Prueba: si, por el contrario, tiene una respuesta lo más precisa posible que no satisface la condición requerida, que significa answer<min(x,y)<=max(x,y)(en cuyo caso min(x,y)es una mejor respuesta, una contradicción) o min(x,y)<=max(x,y)<answer(en cuyo caso max(x,y)es una mejor respuesta, una contradicción).

Así que creo que eso significa que su pregunta se reduce a encontrar la respuesta más precisa posible. Suponiendo aritmética IEEE754 en todo momento, propongo lo siguiente:

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

Mi argumento de que esto da una respuesta más precisa es un análisis de caso algo tedioso. Aquí va:

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

    • Subcase, ni x ni y están desnormalizadas: en este caso, la respuesta calculada x/2.+y/2.manipula las mismas mantisas y, por lo tanto, proporciona exactamente la misma respuesta que el cálculo de (x+y)/2rendiría si asumimos exponentes extendidos para evitar el desbordamiento. Esta respuesta puede depender del modo de redondeo, pero en cualquier caso, IEEE754 garantiza que es la mejor respuesta posible (por el hecho de que lo calculado x+yes una mejor aproximación a x + y matemático, y la división por 2 es exacta en este caso caso).
    • Subcase x está desnormalizado (y así 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.

    • Subcase y está desnormalizado (y así abs(x)>=1): análogo.

  • Caso max(abs(x),abs(y)) < 1.:
    • Subcase, el cálculo x+yno está desnormalizado o está desnormalizado e "par": aunque el cálculo x+ypuede no ser exacto, IEEE754 garantiza que es la mejor aproximación posible a la matemática x + y. En este caso, la división posterior por 2 en la expresión (x+y)/2.es exacta, por lo que la respuesta calculada (x+y)/2.es una mejor aproximación posible a la matemática (x + y) / 2.
    • Subcaso la computado x+yse desnormalizado y "extraño": En este caso exactamente uno de x, y también debe ser desnormalizado-y- "extraño", que significa el otro de X, Y es desnormalizado con el signo opuesto, y por lo que las calculadas x+yes decir exactamente el matemático x + y, (x+y)/2.por lo tanto, IEEE754 garantiza que el cálculo sea la mejor aproximación posible al matemático (x + y) / 2.
Don Hatch
fuente
Me di cuenta de que cuando dije "desnormalizado" realmente quise decir algo más, es decir, números que están tan cerca uno del otro como los números, es decir, el rango de números que es aproximadamente el doble que el rango de números desnormalizados, es decir, los primeros 8 ticks más o menos en el diagrama en en.wikipedia.org/wiki/Denormal_number . El punto es que los "impares" de estos son los únicos números para los cuales dividirlos por dos no es exacto. Necesito reformular esta parte de la respuesta para aclarar esto.
Don Hatch
Fl(opag(X,y))=opag(X,y)(1+δ)El |δEl |tuX/ /2+y/ /2(X+y)/ /2siempre están correctamente redondeados, ausentes sobre / subflujo, todo lo que queda es no mostrar nada sobre- / subflujo, lo cual es fácil.
Kirill
@ Kirill, estoy un poco perdido ... ¿de dónde vienes? Además, no creo que sea del todo cierto que "las divisiones entre 2 son exactas para números no denormales" ... esto es lo mismo que me tropecé, y parece ser un poco incómodo tratar de hacerlo bien. La afirmación precisa es algo más como "x / 2 es exacta siempre que abs (x) sea al menos dos veces el mayor número subnormal" ... ¡argh, incómodo!
Don Hatch
3

Para los formatos de punto flotante binario IEEE-754, ejemplificados por el binary64cálculo (doble precisión), S. Boldo demostró formalmente que el algoritmo simple que se muestra a continuación ofrece el promedio correctamente redondeado.

Sylvie Boldo, "Verificación formal de programas que computan el promedio de coma flotante". En Conferencia Internacional sobre Métodos de Ingeniería Formal , págs. 17-32. Springer, Cham, 2015. ( borrador en línea )

(X+y)/ /2X/ /2+y/ /2binary64C[2-967,2970]C para proporcionar el mejor rendimiento para un caso de uso particular.

Esto produce el siguiente ISO-C99código ejemplar :

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

En un trabajo de seguimiento reciente, S. Boldo y sus coautores mostraron cómo lograr los mejores resultados posibles para los formatos de coma flotante decimal IEEE-754 mediante el uso de operaciones fusionadas de suma múltiple (FMA) y una precisión bien conocida. duplicar el bloque de construcción (TwoSum):

Sylvie Boldo, Florian Faissole y Vincent Tourneur, "Algoritmo formalmente probado para calcular el promedio correcto de números decimales en coma flotante". En el 25 ° Simposio IEEE sobre aritmética informática (ARITH 25) , junio de 2018, pp. 69-75. ( borrador en línea )

njuffa
fuente
2

Aunque puede que no sea un rendimiento súper eficiente, hay una forma muy simple de (1) asegurarse de que ninguno de los números sea mayor que ninguno xo y(sin desbordamientos) y (2) mantener el punto flotante tan "preciso" como posible (y (3) , como un bono adicional, a pesar de que se está usando la resta, no se almacenarán valores como números negativos.

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

De hecho, si realmente desea obtener precisión, ni siquiera necesita realizar la división en el acto; simplemente devuelva los valores de min(x, y)y differenceque puede usar para simplificar lógicamente o manipular más tarde.

IQAndreas
fuente
Lo que estoy tratando de averiguar ahora es cómo hacer que esta misma respuesta funcione con más de dos elementos , manteniendo todas las variables más bajas que el mayor de los números, y usando solo una operación de división para preservar la precisión.
IQAndreas
@becko Sí, estarías haciendo división al menos dos veces. Además, el ejemplo que dio haría que la respuesta saliera mal. Imagine la media de 2,4,9, no es lo mismo que la media de 3,9.
IQAndreas
Tienes razón, mi recursión estaba mal. No estoy seguro de cómo solucionarlo en este momento, sin perder precisión.
Becko
¿Puedes demostrar que esto da el resultado más preciso posible? Es decir, si xy yson punto flotante, su cálculo produce un punto flotante más cercano a (x+y)/2?
Becko
1
¿No se desbordará cuando x, y sean los números expresables mínimo y máximo?
Don Hatch
1

Convierta a mayor precisión, agregue los valores allí y vuelva a convertir.

No debe haber desbordamiento en la precisión más alta y si ambos están en el rango válido de coma flotante, el número calculado también debe estar dentro.

Y debe estar entre ellos, en el peor de los casos, solo la mitad del número mayor si la precisión no es suficiente.

leeroy
fuente
Este es el enfoque de la fuerza bruta. Probablemente funcione, pero estaba buscando un análisis que no requiriera una precisión intermedia más alta. Además, ¿puede estimar cuánta precisión intermedia más alta se requiere? En cualquier caso, no elimine esta respuesta (+1), simplemente no la aceptaré como respuesta.
becko
1

Teóricamente, x/2se puede calcular restando 1 de la mantisa.

Sin embargo, la implementación de operaciones bit a bit como esta no es necesariamente sencilla, especialmente si no conoce el formato de sus números de coma flotante.

Si puede hacer esto, toda la operación se reduce a 3 sumas / restas, lo que debería ser una mejora significativa.

Roland Heath
fuente
0

Estaba pensando en la misma línea que @Roland Heath pero aún no puedo comentar, aquí está mi opinión:

x/2se puede calcular restando 1 del exponente (no la mantisa, restando 1 de la mantisa es restar 2^(value_of_exponent-length_of_mantissa)del valor total).

Sin restricción del caso general, supongamos x < y. (If x > y, vuelva a etiquetar las variables. If x = y, (x+y) / 2es trivial).

  • Transformarse (x+y) / 2en x/2 + y/2, que puede realizarse mediante dos restas de enteros (por una del exponente)
    • Sin embargo, hay un límite inferior en el exponente dependiendo de su representación. Si su exponente ya es mínimo antes de restar 1, este método requerirá un manejo de caso especial. Un exponente mínimo en xse hará x/2más pequeño que representable (suponiendo que mantissa se represente con un 1 implícito).
    • En lugar de restar 1 del exponente de x, desplazar xla mantisa de la derecha por uno (y agregar el primer 1 implícito, si lo hay).
    • Resta 1 del exponente de y, si no es mínimo. Si es mínimo (y es más grande que x, debido a la mantisa), desplace la mantisa a la derecha en uno (agregue el 1 implícito, si lo hay).
    • Desplaza la nueva mantisa de xa la derecha según el exponente de y.
    • Realice la suma de enteros en la mantisa, a menos que la mantisa xse haya desplazado por completo. Si ambos exponentes fueran mínimos, los principales se desbordarán, lo cual está bien, porque se supone que ese desbordamiento se convertirá en uno líder implícito nuevamente.
  • y una adición de coma flotante.
    • No puedo pensar en ningún caso especial aquí; a excepción del redondeo, que también se aplica al desplazamiento descrito anteriormente.
Sin respuesta
fuente