Algoritmo de Kahan: cómo obtener la diferencia exacta de productos

Recientemente volví al análisis de errores de coma flotante para refinar algunos de los detalles en la próxima edición del libro de representación basada en la física . Los números de coma flotante son un área de cálculo interesante, llena de sorpresas (buenas y malas), así como trucos difíciles para deshacerse de las sorpresas desagradables.

En el proceso, me encontré con esta publicación en StackOverflow , de la que aprendí sobre un algoritmo elegante para el cálculo exacto $$$a \ times b-c \ times d$.

Pero antes de continuar con el algoritmo, debe comprender qué es tan astuto en la expresión $$$a \ times b-c \ times d$? Tomar $$$a = 33962.035$, $$$b = -30438.8$, $$$c = 41563.4$y $$$d = -24871.969$. (Estos son los valores reales que obtuve durante el lanzamiento de pbrt ). Para valores flotantes de 32 bits, obtenemos: $$$a \ times b = -1.03376365 \ times 10 ^ 9$y $$$c \ veces d = -1.03376352 \ veces 10 ^ 9$. Llevamos a cabo la resta, y obtenemos $$$-128$. Pero si realiza los cálculos con doble precisión, y al final los convierte en flotantes, obtendrá $$$-75.1656$. Que paso

El problema es que el valor de cada trabajo puede ir mucho más allá del resultado final. $$$-1 \ por 10 ^ 9$, donde la distancia entre valores representables de coma flotante es muy grande - 64. Es decir, al redondear $$$a \ veces b$y $$$c \ veces d$individualmente al flotador representable más cercano, se convierten en números que son múltiplos de 64. A su vez, su diferencia será un múltiplo de 64, y no habrá esperanza de que sea $$$-75.1656$más cerca que $$$-64$. En nuestro caso, el resultado fue aún más debido a cómo se redondearon las dos obras en $$$-1 \ por 10 ^ 9$. Nos encontraremos directamente con la vieja y catastrófica reducción ¹ .

Aquí hay una mejor solución que ² :

inline float DifferenceOfProducts(float a, float b, float c, float d) { float cd = c * d; float err = std::fma(-c, d, cd); float dop = std::fma(a, b, -cd); return dop + err; } 

DifferenceOfProducts() calcula $$$a \ times b-c \ times d$de una manera que evite la contracción catastrófica. Esta técnica fue descrita por primera vez por el legendario William Kahan en el artículo Sobre el costo de la computación de punto flotante sin aritmética extra precisa . Vale la pena señalar que el trabajo de Kahan es interesante de leer en su conjunto, tienen muchos comentarios sobre el estado actual del mundo de los puntos flotantes, así como también consideraciones matemáticas y técnicas. Aquí está uno de sus hallazgos:

Aquellos de nosotros que hemos luchado con las vicisitudes de la aritmética de coma flotante y las "optimizaciones" mal compiladas de los compiladores podemos estar orgullosos de la victoria en esta batalla. Pero si transmitimos la continuación de esta batalla a las generaciones futuras, esto contradecirá toda la esencia de la civilización. Nuestra experiencia dice que los lenguajes de programación y los sistemas de desarrollo son fuentes de demasiado caos con los que tenemos que lidiar. Se pueden prescindir de demasiados errores, así como algunas "optimizaciones" atractivas que son seguras para los enteros, pero que a veces resultan fatales para los números de coma flotante.

Después de rendir homenaje a su ingenio, volvamos a DifferenceOfProducts() : la base del dominio de esta función es el uso de instrucciones FMA ³ fusionadas de suma múltiple. Desde un punto de vista matemático, FMA(a,b,c) es $$$a \ veces b + c$Por lo tanto, al principio parece que esta operación solo es útil como microoptimización: una instrucción en lugar de dos. Sin embargo fma
Tiene una propiedad especial: se completa solo una vez.

En lo habitual $$$a \ veces b + c$primero calculado $$$a \ veces b$, y luego este valor, que en el caso general no se puede representar en formato de coma flotante, se redondea al flotante más cercano. Luego a este valor redondeado se agrega $$$c$, y este resultado se redondea nuevamente al flotador más cercano. FMA se implementa de tal manera que el redondeo se realiza solo al final, un valor intermedio $$$a \ veces b$conserva suficiente precisión, así que después de agregarle $$$c$el resultado final estará más cerca del valor verdadero $$$a \ veces b + c$valor de flotación.

Habiendo tratado con FMA, volveremos a DifferenceOfProducts() . Nuevamente, mostraré las dos primeras líneas:

  float cd = c * d; float err = std::fma(-c, d, cd); 

El primero calcula el valor redondeado $$$c \ veces d$y el segundo ... resta $$$c \ veces d$de su trabajo? Si no sabe cómo funciona FMA, entonces podría pensar que err siempre será cero. Pero cuando se trabaja con FMA, la segunda línea realmente extrae el valor del error de redondeo en el valor calculado $$$c \ veces d$y lo guarda para err . Después de eso, la salida es muy simple:

  float dop = std::fma(a, b, -cd); return dop + err; 

El segundo FMA calcula la diferencia entre los trabajos utilizando el FMA, realizando el redondeo solo al final. En consecuencia, es resistente a la reducción catastrófica, pero necesita trabajar con un valor redondeado. $$$c \ veces d$. La return "parchea" este problema y agrega el error resaltado en la segunda línea. En un artículo de Jeannenrod et al. se muestra que el resultado es verdadero hasta 1.5 ulps (medidas de precisión de la unidad), lo cual es genial: FMA y las operaciones simples de coma flotante son precisas hasta 0.5 ulps, por lo que el algoritmo es casi perfecto.

Usamos un nuevo martillo

Cuando comienzas a buscar formas de usar DifferenceOfProducts() , resulta sorprendentemente útil. Cálculo del discriminante de la ecuación cuadrática? DifferenceOfProducts(b, b, 4 * a, c) llamada de productos DifferenceOfProducts(b, b, 4 * a, c) ⁴ . Cálculo del determinante de una matriz 2x2? El algoritmo resolverá este problema. En la próxima versión de pbrt, encontré alrededor de 80 usos para él. De todos ellos, la función del producto vectorial es la más querida. Siempre ha sido una fuente de problemas, por lo que tuvo que levantar la mano y usar el doble en la implementación para evitar una reducción catastrófica:

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { double v1x = v1.x, v1y = v1.y, v1z = v1.z; double v2x = v2.x, v2y = v2.y, v2z = v2.z; return Vector3f(v1y * v2z - v1z * v2y, v1z * v2x - v1x * v2z, v1x * v2y - v1y * v2x); } 

Y ahora podemos seguir trabajando con float y usar DifferenceOfProducts() .

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { return Vector3f(DifferenceOfProducts(v1.y, v2.z, v1.z, v2.y), DifferenceOfProducts(v1.z, v2.x, v1.x, v2.z), DifferenceOfProducts(v1.x, v2.y, v1.y, v2.x)); } 

Ese astuto ejemplo desde el principio de la publicación es en realidad parte de un trabajo vectorial. En una determinada etapa, el código pbrt necesita calcular el producto vectorial de los vectores $$$(33962.035, 41563.4, 7706.415)$y $$$(- 24871.969, -30438.8, -5643.727)$. Cuando se calcula usando float, obtendríamos un vector $$$(1552, -1248, -128)$. (Regla general: si en los cálculos de coma flotante donde intervienen grandes números, se obtienen valores enteros no tan grandes, entonces esta es una señal casi segura de que se ha producido una reducción catastrófica).

Con doble precisión, el producto vectorial es $$$(1556.0276, -1257.5151, -75.1656)$. Vemos que con flotador $$$x$se ve normal $$$y$ya es bastante malo también $$$z$se convierte en el desastre que se ha convertido en la motivación para encontrar una solución. ¿Y qué resultados obtendremos con DifferenceOfProducts() y valores flotantes? $$$(1556.0276, -1257.5153, -75.1656)$. Valores $$$x$y $$$z$corresponden a doble precisión, y $$$y$ligeramente compensado, de ahí proviene el ulp adicional.

¿Qué hay de la velocidad? DifferenceOfProducts() realiza dos FMA, así como la multiplicación y la suma. Se puede implementar un algoritmo ingenuo con una FMA y una multiplicación, que, al parecer, debería tomar la mitad de tiempo. Pero en la práctica, resulta que después de obtener los valores de los registros, no importa: en el punto de referencia sintético que se encuentra en mi computadora portátil, DifferenceOfProducts() solo 1.09 veces más caro que el algoritmo ingenuo. La operación de doble precisión fue 2.98 veces más lenta.

Tan pronto como se entera de la posibilidad de una reducción catastrófica, todo tipo de expresiones inocentes en el código comienzan a parecer sospechosas. DifferenceOfProducts() parece ser una buena cura para la mayoría de ellos. Es fácil de usar y no tenemos ninguna razón particular para no usarlo.

Notas

1. La reducción catastrófica no es un problema al restar cantidades con diferentes signos o al agregar valores con el mismo signo. Sin embargo, puede convertirse en un problema al agregar valores con diferentes signos. Es decir, las cantidades deben considerarse con la misma sospecha que las diferencias.
2. Como ejercicio para el lector, sugiero escribir la función SumOfProducts() , que proporciona protección contra la contracción catastrófica. Si desea complicar la tarea, explique por qué en DifferenceOfProducts() , dop + err == dop , porque los signos a*b y c*d son diferentes.
3. La instrucción FMA ha estado disponible en la GPU durante más de una década, y en la mayoría de las CPU durante al menos cinco años. Puede ser necesario agregar indicadores del compilador a la CPU para generarlos directamente cuando se usa std::fma() ; en gcc y clang, -march=native trabajos -march=native .
4. En el formato de coma flotante IEEE, la multiplicación por potencias de dos se realiza exactamente, por lo que 4 * a no causa ningún error de redondeo, a menos que se produzca un desbordamiento.