Aproximación rápida de $\tanh(x)$

Question

Esta es una especie de pregunta cruzada de procesamiento de señales/programación/matemáticas. En este momento me parece que está más relacionada con las matemáticas, pero si los moderadores consideran que pertenece a otro lugar, por favor, siéntanse libres de migrarla.

Estoy trabajando en un proyecto en el que tengo una potencia de cálculo limitada y necesito hacer una aproximación rápida de la función tangente hiperbólica sobre un rango bastante grande de argumentos de entrada. Suponiendo que los números se almacenan en punto fijo con una parte fraccionaria de 8 bits entonces la aproximación a $\tanh(x)$ debe trabajar hasta el límite que implica la resolución, o para los argumentos $\tanh^{-1}(\pm[1 - \frac{1}{2^8}]) \approx \pm3.1$ . Sé que las series de Taylor no convergen lo suficientemente rápido en este rango, y no lo he investigado a fondo, pero tampoco creo que una aproximación polinómica de Chebyshev converja lo suficientemente rápido. Por diversas razones tampoco es posible una búsqueda en la tabla.

$\tanh(x)$ es algebraicamente equivalente a $\mathrm{sgn}(x)(1 - \frac{2}{e^{2|x|} + 1})$ . Esto parece prometedor; una expansión en serie de $e^x$ converge mejor que $\tanh(x)$ . Para un cálculo numérico rápido en un hardware limitado, el "+ 1" en el denominador es un obstáculo, ya que significa que hay que realizar una división, y como muchos procesadores sencillos no tienen "división por hardware", la división tiene que hacerse por software mediante repetidos desplazamientos y sustracciones, lo que es terriblemente lento.

Para evitar esto estoy considerando el siguiente enfoque. En lugar de utilizar $\mathrm{sgn}(x)(1 - \frac{2}{e^{2|x|} + 1})$ para representar exactamente tanh(x), defina el siguiente tipo de función similar: $\mathrm{sgn}(x)A(1 - e^{-B|x|})$ para las constantes desconocidas A y B. A continuación, defina una función de pérdida:

$f(x) = \int_0^{3.1}[\tanh(x) - A(1 - e^{-Bx})]^2dx$

A continuación, utilice un método numérico (¿descenso de gradiente?) para encontrar los valores de A y B que minimizan la pérdida. El problema de la aproximación $\tanh(x)$ rápidamente se convierte en una aproximación $e^{-x}$ rápidamente más una resta y una multiplicación extra.

$e^{-x} = \dfrac{1}{e^x} = \dfrac{1}{2^{x \log_2(e)}} = \dfrac{2^{-x}}{2^y}$ donde y es la parte entera y x es la parte fraccionaria. Expresando $e^{-x}$ esta forma me permite reducir el rango del argumento de entrada; dividir un número sin signo por $2^y$ donde y es un número entero positivo, se puede realizar mediante un desplazamiento lógico. Para los argumentos entre 0 y -1 basta con mirar a ojo un Expansión en serie de Taylor de $2^x$ se ve muy bien para sólo 3 términos.

¿Le parece razonable este planteamiento, desde una perspectiva matemática? ¿Problemas que no he tenido en cuenta? Gracias por cualquier consejo.

Edición: Debo añadir que en esta aplicación (procesamiento de audio) estoy dispuesto a sacrificar la precisión absoluta del valor calculado a cambio de la velocidad de procesamiento y de una función que funcione "razonablemente" bien en todo el rango de argumentos de entrada - el valor de retorno para cada argumento en la aproximación no tiene que ser preciso al $2^{-8}$ límite que implica la resolución.

Answer 1

Sugiero el enfoque de Robert Israel. Recientemente he hecho algo similar para aproximar esta función, así que aquí están mis dos centavos. En realidad estaba buscando una aproximación de la función logística, que no es más que una escala y traslación de ésta, pero la tangente hiperbólica es más fácil de trabajar porque es simétrica alrededor del 0. Mi modelo calcula tanh(x/2).

En mi caso encontré una aproximación de grado superior que mantenía el error absoluto máximo por debajo de 2e-7. Me gustaría saber si Maple puede tratar con modelos de orden superior, porque existe el reto de evitar poner polos sobre la línea real. De hecho, este modelo sugerido por Robert tiene un polo aproximadamente en -3.8, que puede no ser un problema en esta aplicación, pero no es bueno tenerlo...

En mis modelos el numerador es un polinomio impar, y el denominador es par, por lo que todos los polos y ceros son pares imaginarios, con un cero extra en el origen. Aquí está el gran modelo que encontré, para un error inferior a 2e-7:

y =  94.9339451088 * x * ((( x**2 + 1.06315869e+03 ) * x**2 + 1.88748783e+05 ) * x**2 + 5.86237309e+06 ) / ((( ( x**2 + 4.03183926e+03 ) * x**2 + 1.64253046e+06 ) * x**2 + 1.28592857e+08 ) * x**2 + 1.11307745e+09 )

Un modelo de menor orden para el error de ~2**-8 también se puede encontrar sin duda con el mismo enfoque. Usé mi código para hacer un nuevo intento de encontrar este modelo más simple, y esto es lo que obtuve, compruébalo:

y = 25.2075466924 * x * ( x**2 + 6.96321678**2 ) / (( x**2 + 3.17572786**2 ) * ( x**2 + 15.57611343**2 ))

Esto aproxima tanh(x/2) con un error ligeramente superior a 2,2e-4 para |x| < 7,6. La derivada de la función aproximada a 7,6 también es igual a 0, por lo que es seguro truncar la función allí (hacerla sgn(x)), mientras que el modelo resultante se convertirá en monotónicamente creciente, una propiedad agradable de tener.

No pude conseguir un buen modelo para un pedido aún más pequeño, pero no creo que pueda ser mucho más simple que eso.

PD: La fracción continuada x / (1+x^2/(3+x^2/(5+...)) también produce una función racional muy similar, pero el error es muy plano alrededor de 0, diferente del ajuste estilo Chebyshev que hice. Me gustaría saber si podríamos utilizar esa fracción continua para llegar a una buena estimación inicial para la regresión no lineal...

Answer 2

Yo escribí este. Es x4 veces más rápido que el construido en tanh(), con sólo 0,0001% de error en todo el rango, y fácil de convertir a las instrucciones vectoriales, como AVX2.

float tanh_c3(float v)
{
    const float c1 = 0.03138777F;
    const float c2 = 0.276281267F;
    const float c_log2f = 1.442695022F;
    v *= c_log2f;
    int intPart = (int)v;
    float x = (v - intPart);
    float xx = x * x;
    float v1 = c_log2f + c2 * xx;
    float v2 = x + xx * c1 * x;
    float v3 = (v2 + v1);
    *((int*)&v3) += intPart << 24;
    float v4 = v2 - v1;
    return (v3 + v4) / (v3 - v4);
}

Answer 3

Métodos seguramente rápidos y precisos para calcular $\ln(x)$ son conocidos. Desde $\tanh^{-1}(x)=\frac12\ln\left(\frac{1+x}{1-x}\right)$ se podría calcular esto en su lugar, y luego ejecutar una búsqueda binaria para encontrar $\tanh(x)$ . Tendrías que ser capaz de calcular con precisión $\ln(Y)$ todo el camino hasta $Y=\frac{1+(1-2^{-8})}{1-(1-2^{-8})}=2^9-1$ . La búsqueda binaria sólo tendría $8$ repeticiones.

Añadido más tarde: Este artículo de Wikipedia describe un método rápido para calcular $\ln(Y)$ utilizando la media aritmética-geométrica.

Answer 4

Como tú dices,

$$\tanh(x)=\operatorname{sgn}(x)\left(1-\frac2{e^{2|x|}+1}\right)$$

Por expansión geométrica podemos reescribir esto como

$$\tanh(x)=\operatorname{sgn}(x)\left(1-2e^{-2|x|}+2e^{-4|x|}-2e^{-6|x|}+\dots\right)$$

Esta expansión converge rápidamente con un error de orden $e^{-2N|x|}$ para $N$ términos, que es pequeño cuando $x$ es grande. Evita cualquier división aquí, y todos los exponenciales se reducen a uno ya que $e^{-2n|x|}=(e^{-2|x|})^n$ .

A continuación mostramos un gráfico de las aproximaciones hasta el $\pm2e^{-2n|x|}$ plazo hasta $n=8$ :

Lo único que queda por hacer es aproximarlo cerca del origen, lo que puede hacerse fácilmente mediante aproximaciones de Taylor.

Si esto sigue siendo insatisfactorio, se puede desear una aproximación alternativa para $x$ entre, por ejemplo, 0,5 y 1.

Answer 5

¿Has mirado CORDIC ? (por desgracia, no he visto el sin búsqueda de tablas pero lo dejaré aquí de todos modos)