Cancelación catastrófica en logsum

18

Estoy tratando de implementar la siguiente función en coma flotante de doble precisión con un error relativo bajo :

logsum(x,y)=log(exp(x)+exp(y))

Esto se usa ampliamente en aplicaciones estadísticas para agregar probabilidades o densidades de probabilidad que se representan en el espacio logarítmico. Por supuesto, exp(x) o exp(y) podrían desbordarse o desbordarse fácilmente, lo que sería malo porque el espacio de registro se usa para evitar el desbordamiento en primer lugar. Esta es la solución típica:

logsum(x,y)=x+log1p(exp(yx))

La cancelación de yx ocurre, pero se mitiga con exp . Peor con diferencia es cuando x y log1p(exp(yx)) están cerca. Aquí hay un gráfico de error relativo:

ingrese la descripción de la imagen aquí

El gráfico se corta entre para enfatizar la forma de la curva l o g s u m ( x , y ) = 0 , alrededor de la cual ocurre la cancelación. He visto el error de hasta 10 - 11 y sospechan que se pone mucho peor. (FWIW, la función de "verdad sobre el terreno" se implementa utilizando los flotantes de precisión arbitraria de MPFR con una precisión de 128 bits).1014logsum(x,y)=01011

He intentado otras reformulaciones, todas con el mismo resultado. Con como la expresión externa, se produce el mismo error al tomar un registro de algo cercano a 1. Con l o g 1 p como la expresión externa, la cancelación ocurre en la expresión interna.loglog1p

Ahora, el error absoluto es muy pequeño, por lo que tiene un error relativo muy pequeño (dentro de un épsilon). Uno podría argumentar que, debido a que un usuario de l o g s u m está realmente interesado en las probabilidades (no en las probabilidades de registro), este terrible error relativo no es un problema. Es probable que generalmente no lo sea, pero estoy escribiendo una función de biblioteca, y me gustaría que sus clientes puedan contar con un error relativo no mucho peor que el error de redondeo.exp(logsum(x,y))logsum

Parece que necesito un nuevo enfoque. ¿Qué puede ser?

Neil Toronto
fuente
No entiendo tu último párrafo. "dentro de un épsilon" no significa nada para mí. ¿Te refieres a una Unidad en el Último Lugar ? En cuanto a los usuarios interesados ​​en las probabilidades, un error de probabilidad de registro pequeño dará como resultado un error de probabilidad grande, por lo que este no es el caso.
Aron Ahmadia
Por curiosidad, ¿has intentado tomar el "mejor" de tus dos métodos y trazar el error de eso? Entonces, todo lo que necesita es la lógica correcta para detectar en qué caso se encuentra (es de esperar que sea menos costoso o parte del costo requerido del algoritmo de todos modos), luego cambie al método apropiado.
Aron Ahmadia
@AronAhmadia: "Dentro de un épsilon" significa un error relativo menor que un épsilon de coma flotante de doble precisión, que es aproximadamente 2.22e-16. Para flotadores normales (es decir, no subnormales), corresponde a aproximadamente una ulp. Además, si es el error absoluto de x , entonces el error relativo de exp ( x ) es exp ( a ) - 1 , que es casi la función de identidad cerca de cero. IOW, un pequeño error absoluto para x implica un pequeño error relativo para exp (axexp(x)exp(a)1x . exp(x)
Neil Toronto
Anexo: Cuando el error absoluto está cerca de cero. Cuando a > 1 , por ejemplo, tienes razón: el relativo explota. aa>1
Neil Toronto

Respuestas:

12

La fórmula debe ser numéricamente estable. Se generaliza a un cálculo numéricamente estable de Iniciar sesión

logsum(x,y)=max(x,y)+log1p(exp(abs(xy))
logiexi=ξ+logiexiξ,   ξ=maxixi

logsum(x,y)=max(x,y)+lexp(xy)
lexp(z):=log(1+e|z|)
z
Arnold Neumaier
fuente
En términos de error absoluto, lo es. En términos de error relativo, es horrible cuando la salida está cerca de cero.
Neil Toronto
xy
Para x = -0.775 e y = -0.6175, obtengo un error de 62271 ulps y un error relativo de 1.007e-11.
Neil Toronto
1
Calcule puntos de datos altamente precisos en el rango de interés: se necesitan al menos dos rangos diferentes debido al comportamiento asintótico. Se puede usar la expresión definitoria para z no cercana a cero. Para el rango excepcional, ajuste una función racional de grado suficientemente alto para obtener la precisión deseada. Para la estabilidad numérica, use polinomios de Bernstein o polinomios de Tchebychev en numerador y denominador, adaptados al intervalo de interés. Al final, expanda a una fracción continua y descubra cuánto se pueden truncar los coeficientes sin afectar la precisión.
Arnold Neumaier
1
l=l(z)m