Interpolación lineal precisa de punto flotante

Question

Quiero realizar una interpolación lineal simple entre $A$ y $B$ (que son valores binarios de punto flotante) utilizando matemáticas de punto flotante con reglas de redondeo IEEE-754 al más cercano o al más cercano par, con la mayor precisión posible. Tenga en cuenta que la velocidad no es una gran preocupación aquí.

Conozco dos enfoques básicos. Usaré los símbolos $\oplus, \ominus, \otimes, \oslash$ siguiendo a Knuth [1], para significar suma, resta, producto y división de punto flotante, respectivamente (en realidad no uso la división, pero la he enumerado por completitud).

(1) $\quad f(t) = A\,\oplus\,(B\ominus A)\otimes t$

(2) $\quad f(t) = A\otimes(1\ominus t)\,\oplus \,B\otimes t$

Cada método tiene sus pros y sus contras. El método (1) es claramente monótono, lo cual es una propiedad muy interesante, mientras que no me resulta obvio en absoluto que eso sea cierto para el método (2), y sospecho que puede que no sea el caso. Por otro lado, el método (2) tiene la ventaja de que cuando $t = 1$ el resultado es exactamente $B$, no una aproximación, y esa es también una propiedad deseable (y exactamente $A$ cuando $t = 0$, pero el método (1) también hace eso). Eso se sigue de las propiedades enumeradas en [2], en particular:

$u\oplus v = v\oplus u$

$u\ominus v = u\oplus -v$

$u\oplus v = 0$ si y solo si $v = -u$

$u\oplus 0 = u$

$u\otimes 1 = u$

$u\otimes v = 0$ si y solo si $u = 0$ o $v = 0$

En [3] Knuth también discute este caso:

$u' = (u\oplus v)\ominus v$

lo que implica que $u'$ puede o no ser igual a $u$. Reemplazando $u$ con $B$ y $v$ con $-A$ y usando las reglas anteriores, se sigue que no está garantizado que $A\oplus(B\ominus A) = B$, lo que significa que el método (1) no siempre produce $B$ cuando $t = 1$.

Entonces, aquí vienen mis preguntas:

1. ¿Está garantizado que el método (2) es monótono?
2. Si no, ¿hay un método mejor que sea preciso, monótono y produzca $A$ cuando $t = 0$ y $B$ cuando $t = 1$?
3. Si no (o no lo sabe), ¿el método (1) cuando $t = 1$ siempre se excede (es decir, $A\oplus(B\ominus A)=A+(B-A)\cdot t$ para algún $t \geq 1$)? ¿Siempre se queda corto (idem para algún $t \leq 1$)? ¿O a veces se excede y a veces no llega?

Supongo que si el método (1) siempre se queda corto, puedo hacer un caso especial cuando $t = 1$ para obtener la propiedad deseada de ser exactamente igual a $B$ cuando $t = 1$, pero si siempre se excede, entonces no puedo. Esa es la razón de la pregunta 3.

EDITAR: He encontrado que la respuesta a la pregunta 3 es que a veces se excede y a veces no llega. Por ejemplo, en doble precisión:

-0x1.cae164da859c9p-1 + (0x1.eb4bf7b6b2d6ep-1 - (-0x1.cae164da859c9p-1))

da como resultado 0x1.eb4bf7b6b2d6fp-1, que es 1 ulp mayor que el original, mientras que

-0x1.be03888ad585cp-1 + (0x1.0d9940702d541p-1 - (-0x1.be03888ad585cp-1))

da como resultado 0x1.0d9940702d540p-1, que es 1 ulp menor que el original. Por otro lado, el método que planeé (caso especial de $t=1$) no funcionará, porque he encontrado que puede ser el caso donde $t < 1$ y $A\oplus(B\ominus A)\otimes t > B$, por ejemplo:

~~

t = 0x1.fffffffffffffp-1
A = 0x1.afb669777cbfdp+2
B = 0x1.bd7b786d2fd28p+1

$A \oplus (B \ominus A)\otimes t =\,$ 0x1.bd7b786d2fd29p+1

~~

lo que significa que si se va a utilizar el método (1), la única estrategia que puede funcionar es recortar.

Actualización: Como señaló Davis Herring en un comentario y verificado posteriormente por mí, hacer un caso especial de t=1 funciona en realidad.

---

Referencias

[1] D.E.Knuth, The Art of Computer Programming, vol. 2: Seminumerical algorithms, tercera edición, pág. 215

[2] Ibíd. pp. 230-231

[3] Ibíd. p.235 ec.(41)

Answer 1

1) El método 2 no siempre es monótono. Contraejemplo en doble precisión:

$$A = B = 4000\quad t=\tt{0x1.b2f3db7800a39p-2}$$

A partir de ahí, resulta que:

$$1\ominus t=\tt{0x1.26861243ffae4p-1}$$

$$4000\otimes \tt{0x1.26861243ffae4p-1}\oplus 4000\otimes \tt{0x1.b2f3db7800a39p-2}=\tt{0x1.f400000000001p+11}\approx 4000.0000000000005$$

(obviamente, cuando $t = 0$ y $t = 1$ el resultado es igual a 4000, por lo que asciende y luego desciende, por lo tanto no es monótono).

2) La segunda pregunta pide un mejor método. Aquí hay un candidato:

$$lerp_{A,B}(t)=\begin{cases}A\oplus(B\ominus A)\otimes t&& \text{if}\,\ t<0.5 \\ B\ominus(B\ominus A)\otimes(1\ominus t)&&\text{otherwise}\end{cases}$$

Este método tiene las siguientes propiedades:

• Coincide con los extremos (obvio).
• Ambas mitades son monótonas (obvio).
• Exacto cuando $A=B$ (obvio)
• Parece ser monótono en el cambio de intervalos.

Este último es la única parte que no está probada. Resulta que hay muchos valores de $A$ y $B$ con $AB\ominus(B\ominus A)\otimes 0.5$; sin embargo, después de probar muchos miles de millones de números en precisión simple y doble, no pude encontrar un solo caso que incumpla ninguna de estas implicaciones:

Sea $s=0.5-\mathrm{ulp}(0.5)/2\quad$(el número inmediatamente precedente a 0.5), entonces

$A

$A>B\implies A\oplus(B\ominus A)\otimes s\ge B\ominus(B\ominus A)\otimes 0.5$

Actualización: La razón por la que funciona parece tener que ver con este hecho: excluyendo el desbordamiento, para cualquier número binario de punto flotante positivo $u, u\otimes s < u/2$, y de manera similar para $u$ negativo cambiando la dirección de la desigualdad.

O dicho de otra manera, bajo la misma suposición, para base 2 y precisión $p, u\otimes\frac{(2^p-1)}{2^p}

---

He eliminado mis otras dos propuestas después de descubrir que ambas violaban la monotonía.

Answer 2

El valor interpolado exacto puede ser representado como la suma aritmética de 5 valores de punto flotante: uno para el resultado aproximado, uno para el error en cada operación [1] y uno para el error del producto al multiplicar el error de sustracción por t. Aquí hay una versión en C del algoritmo en precisión simple; para simplificar, utiliza la multiplicación-adición fusionada para calcular el error del producto:

#include 

float exact_sum(float a, float b, float *err)
{
    float sum = a + b;
    float z = sum - a;
    *err = a - (sum - z) + (b - z);
    return sum;
}

float exact_mul(float a, float b, float *err)
{
    float prod = a * b;
    *err = fmaf(a, b, -prod);
    return prod;
}

float exact_lerp(float A, float B, float t,
                 float *err1, float *err2, float *err3, float *err4)
{
    float diff = exact_sum(B, -A, err1);
    float prod = exact_mul(diff, t, err2);
    *err1 = exact_mul(*err1, t, err4);
    return exact_sum(A, prod, err3);
}

Para que este algoritmo funcione, las operaciones deben cumplir con las normas IEEE-754 en el modo de redondeo al más cercano o par. Eso no está garantizado por el estándar de C, pero el compilador GNU gcc puede ser instruido para hacerlo, al menos en procesadores que admiten SSE2 [2][3].

Se garantiza que la adición aritmética (a diferencia de la adición de punto flotante) de $resultado + err_1 + err_2 + err_3 + err_4$ será igual al resultado deseado; sin embargo, no se garantiza que incluso sea un número de punto flotante. No estoy al tanto de ninguna garantía para las sumas parciales de estos. Se espera que $err_4$ sea generalmente mucho más pequeño que los demás cuando no son cero.

Por ejemplo: exact_lerp(0.23456789553165435791015625f, 7.345678806304931640625f, 0.300000011920928955078125f, &err1, &err2, &err3, &err4) devuelve $2.3679010868072509765625$ y los errores son $6.7055225372314453125\cdot 10^{-8}$, $8.4771045294473879039287567138671875\cdot 10^{-8}$, $1.490116119384765625\cdot 10^{-8}$ y $2.66453525910037569701671600341796875\cdot 10^{-15}$ respectivamente. Estos números suman el resultado exacto, que es $2.36790125353468550173374751466326415538787841796875$ (no es un flotante de precisión simple, aunque en este caso cabe en uno de doble precisión).

Todos los números en el ejemplo anterior están escritos con sus valores exactos, en lugar de un número que se aproxima a ellos. Todos menos el último son flotantes de precisión simple.

Si el objetivo es precisión, esperaría que calcular $err_1 \oplus err_2 \oplus err_3 \oplus err_4 \oplus resultado$ (en orden de izquierda a derecha) dé una mejor aproximación al resultado real que solo usar $resultado$ por sí solo. En el ejemplo anterior, da $2.367901325225830078125$ que coincide con el valor redondeado del resultado exacto. No he realizado ninguna prueba en cuanto a la monotonicidad. $t=0$ coincide con el punto final, y probablemente $t=1$ también lo hace.

No estoy seguro de si vale la pena agregar $err_4$ en absoluto, dado lo pequeña que es su contribución.

Referencias

• [1] Graillat, Stef (2007). Producto y Exponentiación Precisos en Punto Flotante.
• [2] https://stackoverflow.com/questions/7295861/enabling-strict-floating-point-mode-in-gcc
• [3] Semántica de la Matemática en Punto Flotante en GCC

Answer 3

Mostraré cómo calcular el valor de la interpolación lineal con redondeo correcto; es decir, como si la interpolación lineal se calculara con aritmética exacta y luego se redondeara.

Primero represente la interpolación lineal como $a + tb - ta$. Utilizando el producto exacto de la respuesta de Pedro Gimenos, descomponga $tb = c_0 + c_1$. De manera similar para $ta = d_0 + d_1$. Estas transformaciones no tienen errores.

El problema se reduce entonces a calcular la suma $a + c_0 + c_1 + d_0 + d_1$ con redondeo correcto. Esto se puede hacer aplicando el algoritmo iFastSum del artículo Correct Rounding and a Hybrid Approach to Exact Floating-Point Summation.

Como bonificación, el redondeo correcto implica que el signo del resultado es exacto.

Claramente, el mismo procedimiento se puede adaptar para calcular productos punto con redondeo correcto. Esto implica que hay una prueba exacta y eficiente del signo de un producto punto, lo cual es bastante interesante para la geometría computacional.

Hice una rápida medición de tiempo. El algoritmo descrito aquí, implementado en C#, es 14 veces más lento que el ingenuo punto flotante $(1 - t)a + tb$, promediado en 10 millones de cálculos donde $a$, $b$ y $t$ son cada uno uniformemente aleatorios en $[0, 1]$.

Nota: para un algoritmo de sumado correctamente redondeado que es rápidamente asintótico y utiliza iFastSum como subalgoritmo: "Algorithm 908: Online Exact Summation of Floating-Point Streams", ACM Transactions on Mathematical Software,Volume 37, Issue 3, September 2010.

Answer 4

El Método 2 es monótono, si estás usando redondeo al más cercano hacia el par (que afortunadamente es el valor predeterminado).

Vamos a considerar $A=B=1$ y números de media precisión (porque son cortos), y $t=5/2^{12}$: t = 0.000000000101 1-t = 0.111111111011 Pero eso es demasiado largo - solo obtenemos 11 bits. El valor que realmente obtenemos depende del modo de redondeo en el que estemos.

en los modos "hacia 0" (truncar) y "hacia $-\infty$" (piso):

1-t = 0.11111111101

en los otros modos, "hacia $\infty$" (techo), "redondeo lejos de 0", y "al par":

1-t = 0.11111111110

Ahora vamos a sumarlos de nuevo.

truncar y piso: $t+(1-t)=0.11111111111(1)=0.11111111111<1$

techo y redondeo lejos de 0: $t+(1-t)=1.0000000000(1) = 1.0000000001>1$

al par: $t+(1-t)=1.0000000000(1)=1.0000000000=1$

Un análisis adicional nos dice qué está pasando: el objetivo es que los dos pasos de redondeo se contrarresten mutuamente. Esto nunca sucede con piso/truncar/techo. La mayor parte del tiempo sucede con redondeo lejos de cero, pero en la situación donde hay un empate, ambos pasos de redondeo sesgan el resultado hacia arriba. Con redondeo al más cercano, los pasos de redondeo siempre son opuestos entre sí: para aquellos que redondean hacia abajo durante el paso de $1-t$ (por ejemplo, $3/2^{12}$), redondearán hacia arriba durante el paso de la adición, y viceversa.

Answer 5

En mi otra respuesta proporcioné un enfoque correctamente redondeado. Aquí proporcionaré otro enfoque para una precisión mejorada que está dentro del 2x rendimiento de la interpolación ingenua, y probablemente te dará un resultado fielmente redondeado (es decir, error < 1 ulp).

Simplemente implementa la aritmética doble de pares compensados del artículo "Cálculos de Punto Flotante Redondeados Fielmente". Luego realiza el cálculo como lo harías normalmente. El esfuerzo de implementación es bajo.