Es bien sabido que cuando x está cerca de 0, el cálculo de exp(x) - 1 introduce importantes errores de cancelación. Por ello, tenemos expm1 implementado en c99 y python.
Mi pregunta es cómo expm1 evita el error de cancelación. ¿Puede alguien darme una idea general sin demasiados detalles de implementación?
La razón de hacer la función expm1 es evitar la pérdida de precisión al calcular $e^x-1$ cuando $x$ está cerca de $0$ . Para los grandes $x$ o muy negativo $x$ la función expm1 dará el mismo resultado que calcular la exp habitual y luego restar $1$ por lo que sólo nos centraremos en $x\approx 0$ donde se encuentran las diferencias. Pero puedo mencionar que una forma común de calcular exp (y algo similar funciona para expm1) para grandes $|x|$ es escribir $e^{x} = (e^{x/2})^2 = ((e^{x/4})^2)^2$ y así sucesivamente (este truco se explica más en la otra respuesta). Esto nos permite calcular exp para un argumento mucho más pequeño, para el que funciona bien una aproximación tipo serie de Taylor, y luego hacer unas cuantas multiplicaciones para recuperar el valor que queremos.
Ahora para explicar lo que sucede cerca de $x = 0$ utilicemos la forma más sencilla de calcular $e^x$ para los más pequeños $x$ , a saber ${\rm exp}(x) = 1 + x$ como ejemplo. Asimismo, nuestra implementación de expm1 será ${\rm expm1}(x) = x$ . Esto no está muy lejos de cómo lo haría una aplicación real para los más pequeños $x$ .
Para entender por qué expm1 $(x)$ es superior a exp $(x)-1$ cerca de $x = 0$ Vamos a repasar rápidamente cómo funcionan los números de punto flotante. Tienen la forma $a\cdot 2^b$ donde tenemos un cierto número de bits para representar $a$ y $b$ (aquí $a$ es un número entre $1$ y $2$ y $b$ algún número entre $-N$ y $N$ ). Para el caso de los números de doble precisión $a$ tiene una precisión de $10^{-16}$ y $N\sim 1023$ . Esto es muy simplificado, pero es todo lo que necesitamos aquí.
El número más pequeño que podemos representar arriba $1$ será $\approx 1 + 10^{-16}$ (el exponente $b=0$ por lo que la precisión es la precisión de $a$ ). Por otro lado el número más pequeño próximo a cero que podemos representar es $2^{-1023} \approx 10^{-308}$ (aquí $a=1$ y la precisión viene determinada por $b$ ). Este es el punto clave a entender: los números que podemos representar no están espaciados uniformemente (como $0,\epsilon,2\epsilon,3\epsilon,\ldots, 1-\epsilon,1,1+\epsilon, \ldots$ ) pero el tamaño del espacio entre dos números que podemos representar depende del número que estemos mirando y es mucho más pequeño cerca de $0$ que cerca de $1$ .
Por ello, si $x$ es menor que $10^{-16}$ entonces $1+x$ simplemente se evalúa como $1$ y $(1+x) - 1$ evalúa a $0$ . Si queremos calcular $\frac{e^x-1}{x}$ utilizando el $e^x$ por lo tanto, esto dará un resultado muy erróneo $\frac{0}{x} = 0$ . Si por el contrario utilizamos la función expm1 entonces tenemos $\frac{(e^x-1)}{x} = \frac{x}{x} = 1$ lo cual es correcto (por supuesto si $|x|<10^{-308}$ entonces esto también dará $0$ ya que no podemos representar números más pequeños).
Esto era sólo un ejemplo sencillo, pero explica el problema principal: obtenemos enormes errores de cancelación al añadir $1$ a algo pequeño y luego restar $1$ al calcular $e^x-1$ . Para saber más sobre cómo se aplica realmente, consulte esta nota y el Implementación de Google Go (estos enlaces fueron dados en los comentarios de David K y Bungo).
En $x=0$ el exponencial se calcula como una versión de $1+xp(x)$ o $\frac{1+xg(x^2)}{1-xg(x^2)}$ o similar con algunos polinomios que aproximan el término exacto detrás de ellos de manera uniforme en algún intervalo.
Para otros valores las leyes del logaritmo $e^x=(e^{x/2^m})^{2^m}$ o $e^x=2^k3^me^{x-k\ln2-m\ln3}$ se utilizan para reducir el argumento a algún valor pequeño.
Entonces es fácil modificar esas expresiones de aproximación para que devuelvan el expm1(x) valor, en los casos mencionados $xp(x)$ o $\frac{2xg(x^2)}{1-xg(x^2)}$ .
I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.
