Comparación relativa de números de coma flotante

10

Tengo una función numérica que f(x, y)devuelve un número de coma flotante doble que implementa alguna fórmula y quiero verificar que sea correcta contra las expresiones analíticas para todas las combinaciones de los parámetros xy yque me interese. ¿Cuál es la forma correcta de comparar lo calculado y números analíticos de coma flotante?

Digamos que los dos números son ay b. Hasta ahora me he asegurado de que tanto los errores absolutos ( abs(a-b) < eps) como los relativos ( abs(a-b)/max(abs(a), abs(b)) < eps) sean menores que eps. De esa forma, detectará imprecisiones numéricas incluso si los números son, digamos, alrededor de 1e-20.

Sin embargo, hoy descubrí un problema, el valor numérico ay el valor analítico bfueron:

In [47]: a                                                                     
Out[47]: 5.9781943146790832e-322

In [48]: b                                                                     
Out[48]: 6.0276008792632078e-322

In [50]: abs(a-b)                                                              
Out[50]: 4.9406564584124654e-324

In [52]: abs(a-b) / max(a, b)                                                  
Out[52]: 0.0081967213114754103

Entonces, el error absoluto [50] es (obviamente) pequeño, pero el error relativo [52] es grande. Entonces pensé que tenía un error en mi programa. Al depurar, me di cuenta, que estos números son normales . Como tal, escribí la siguiente rutina para hacer la comparación relativa adecuada:

real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
    r = 0
else
    r = d / m
end if
end function

Donde tiny(1._dp)devuelve 2.22507385850720138E-308 en mi computadora. Ahora todo funciona y simplemente obtengo 0 como error relativo y todo está bien. En particular, el error relativo anterior [52] es incorrecto, simplemente es causado por una precisión insuficiente de los números denormales. ¿Es rel_errorcorrecta mi implementación de la función? ¿Debo comprobar que abs(a-b)es menor que minúsculo (= denormal) y devolver 0? ¿O debería comprobar alguna otra combinación, como max(abs(a), abs(b))?

Me gustaría saber cuál es la forma "adecuada".

Ondřej Čertík
fuente

Respuestas:

11

Puede comprobar directamente para denormals utilizando isnormal()partir de math.h(C99 o temprano, POSIX.1 o posterior). En Fortran, si el módulo ieee_arithmeticestá disponible, puede usarlo ieee_is_normal(). Para ser más precisos acerca de la igualdad difusa, debe considerar la representación en coma flotante de denormales y decidir qué quiere decir para que los resultados sean lo suficientemente buenos.

Más aún, para creer que cualquiera de los resultados es correcto, debe asegurarse de no perder demasiados dígitos en un paso intermedio. La computación con denormales generalmente no es confiable y debe evitarse haciendo que su algoritmo se vuelva a escalar internamente. Para asegurarse de que su escala interna fue exitosa, le recomiendo activar las excepciones de punto flotante usando feenableexcept()C99 o el ieee_arithmeticmódulo en Fortran.

Aunque puede hacer que su aplicación capte la señal que se genera en excepciones de coma flotante, todos los núcleos que he intentado restablecen el indicador de hardware para fetestexcept()que no devuelva un resultado útil. Cuando se ejecuta con -fp_trap, los programas PETSc (por defecto) imprimirán un seguimiento de la pila cuando se genera un error de coma flotante, pero no identificarán la línea ofensiva. Si se ejecuta en un depurador, el depurador conserva el indicador de hardware y rompe la expresión ofensiva. Puede verificar el motivo exacto llamando fetestexceptdesde el depurador donde el resultado es un OR bit a bit de los siguientes indicadores (los valores pueden variar según la máquina, vea fenv.h; estos valores son para x86-64 con glibc).

  • FE_INVALID = 0x1
  • FE_DIVBYZERO = 0x4
  • FE_OVERFLOW = 0x8
  • FE_UNDERFLOW = 0x10
  • FE_INEXACT = 0x20
Jed Brown
fuente
Gracias por la excelente respuesta. La expresión analítica con la que comparo en el régimen asintótico es exp(log_gamma(m+0.5_dp) - (m+0.5_dp)*log(t)) / 2para m = 234, t = 2000. Va a cero rápidamente a medida que aumento m. Todo lo que quiero es asegurarme de que mi rutina numérica devuelva números "correctos" (devolver cero también está perfectamente bien) a al menos 12 dígitos significativos. Entonces, si el cálculo devuelve un número normal, entonces es simplemente cero y no debería haber ningún problema. Entonces, solo la rutina de comparación debe ser robusta contra esto.
Ondřej Čertík
5

Donald Knuth tiene una propuesta para un algoritmo de comparación de coma flotante en el volumen 2 "Algoritmos seminéricos" de "El arte de la programación de computadoras". Fue implementado en C por Th. Belding (ver paquete fcmp ) y está disponible en GSL .

GertVdE
fuente
2
Aquí está mi implementación de Fortran: gist.github.com/3776847 , tenga en cuenta que de todos modos necesito manejar explícitamente números denormales. De lo contrario, creo que es bastante equivalente al error relativo, la única diferencia es que, en lugar de hacerlo abs(a-b)/max(a, b) < eps, lo hacemos abs(a-b)/2**exponent(max(a, b)) < eps, lo que prácticamente deja caer la mantisa en el max(a, b), por lo que, en mi opinión, la diferencia es insignificante.
Ondřej Čertík
5

Los números desnormalizados óptimamente redondeados pueden tener un alto error relativo. (Vaciar eso a cero mientras todavía lo llama un error relativo es engañoso).

Pero cerca de cero, calcular los errores relativos no tiene sentido.

Por lo tanto, incluso antes de alcanzar números desnormalizados, probablemente debería cambiar a una precisión absoluta (es decir, la que desea garantizar en este caso).

yXEl |y-XEl |unasisunaCC+rmilunaCCmax(El |XEl |,El |yEl |)

Entonces, los usuarios de su código saben con precisión cuánta precisión tienen realmente.

Arnold Neumaier
fuente
¿Estás seguro de que no tiene sentido calcular errores relativos cercanos a cero? Creo que no tiene sentido solo si hay pérdida de precisión (por cualquier razón). Si, por ejemplo, hay pérdida de precisión para x <1e-150 debido a algunos problemas numéricos (como restar dos números grandes), entonces tiene razón. En mi caso, sin embargo, los números parecen ser precisos hasta cero, excepto cuando llega a los números denormales. Entonces, en mi caso, absacc = 1e-320 más o menos y puedo verificar abs(a-b) < tiny(1._dp)como lo hice anteriormente.
Ondřej Čertík
@ OndřejČertík: en ese caso, reemplace el 1e-150 por 1e-300 o cualquier límite que pueda verificar. En cualquier caso, muy cerca de cero, comete un error absoluto, y su reclamo de error debe reflejar esto en lugar de declarar que el error relativo es cero.
Arnold Neumaier
Veo. Puedo verificar que todo funciona para números superiores a tiny(1._dp)=2.22507385850720138E-308(cometí un error en mi comentario anterior, es 2e-308, no 1e-320). Entonces este es mi error absoluto. Entonces necesito comparar el error relativo. Veo tu punto, creo que tienes razón. ¡Gracias!
Ondřej Čertík
1
El |y-XEl |-unasisunaCCmax(El |XEl |,El |yEl |)