Encontrar una buena aproximación para$x^{1/2.4}$

Question

Me gustaría una buena (8 bits de precisión) la aproximación de $x^{1/2.4}$ en el rango $[0, 1]$. Esta transformación se utiliza para la conversión lineal de las intensidades a SRGB comprimido valores, por lo que es importante que puedo hacer es correr rápido.

La trama de la función:

El uso de un simple polinomio no es práctico porque

1. la función tiene un montón de alta derivadas de orden
2. la función es de aproximadamente asintótica $x$, que es muy diferente el comportamiento asintótico de los polinomios de alto orden

Yo ya tengo el código que construye un arbitrario de grado del polinomio para cualquier función de minimizar el error de los cuadrados e incluso para un 10mo grado del polinomio, la exactitud es todavía sólo como 6 bits significativos.

Luego me enteré de la función racional aproximaciones, lo que tendrá un mejor comportamiento asintótico. Pero el problema es que no sé cómo encontrar el óptimo de los coeficientes. Hay el Pade formulación que crea una aproximación en torno a un solo punto, pero ya que no uso global de la información, puede ser un muy mal ajuste en general como en series de Taylor.

Yo había Mathematica crear una aproximación de la forma $(a_0 + a_1 x + a_2 x^2) / (1 + a_3 x)$ con PadeApproximant[x^(1/2.4), {x, 0.2, {3, 2}}], que es mucho mejor que un simple 3er grado del polinomio, pero todavía no lo suficientemente bueno, así que quiero encontrar a un nivel global de la solución óptima, probablemente de la misma forma.

Traté de encontrar una solución de mínimos cuadrados como antes, pero implica 4 enormes, no-lineal de ecuaciones polinómicas, que se está llevando a Mathematica para siempre (he esperado 1/2 hora hasta el momento) para resolver.

Puede alguien sugerir cómo resolver esas ecuaciones no lineales, o de otra manera de encontrar una función racional de aproximación, o de una forma totalmente diferente de aproximación?

Gracias por la ayuda

Answer 1

Yo creo que tu razonamiento es correcto, una aproximación racional es probablemente la mejor solución. Pero, la construcción de aproximaciones racionales es difícil.

Usted tiene que decidir cómo va a medir el error. Si usted acaba de ver en la máxima diferencia entre el valor de la función y de la aproximación, entonces usted está haciendo "uniforme" o "minimax" aproximación, que por lo general es más difícil que el de mínimos cuadrados aproximación. Como usted ha mencionado, Pade aproximación no es muy bueno para esto, porque es una Taylor-serie-como en el enfoque que se fija en un solo punto.

Yo recomendaría que pruebes el Mathematica RationalInterpolation función, si no lo has hecho ya.

O, si usted tiene acceso a Matlab, hay un add-on llamado Chebfun que hace un muy buen trabajo de construcción de minimax polinómicas y racionales aproximaciones. Hay comandos llamado ratinterp y remez, y un par de otros. El nombre de "remez" viene de la Remez Cambio de Algoritmo, que es la forma estándar de la informática minimax aproximaciones. Esta sección de la Chebfun documentación explica las técnicas que están disponibles. Sección 4.8 cubre aproximaciones racionales.

Hay un buen libro por Trefethen que ofrece una moderna cómputo orientado a la cuenta de la teoría de la aproximación, con un enfoque en Chebfun. Intuitivo y fácil de leer.

Editar

He probado el siguiente código de Mathematica

gx = GeneralMiniMaxApproximation[{t, 1 - t + t^(1/2.4)}, 
                             {t, {0, 1}, 3, 3}, x, 
                              MaxIterations -> 200, WorkingPrecision -> 50]

y, a continuación,

hx = gx - 1 + x

Esto le dio una aproximación a $hx$ de su función

$$ \frac{x^4 - 1.3730926529 x^3 - 1.2231795079 x^2 - 0.0120332034 x - 3.892203211 \times 10^{-7}} {x^3 - 2.47604929086 x^2 - 0.13996763348 x - 0.00008072719}$$

La función de error se parece a esto:

No se muy bien 8 bits, pero muy cerca de llegar. Usted puede jugar con los grados del numerador y el denominador para intentar conseguir algo mejor.

Answer 2

Considere la posibilidad de $\frac{1}{2.4} \approx 0.41$, tal vez comenzando con una raíz cuadrada y la aproximación de la proporción es más fácil. O simplemente el uso de aproximaciones para el logaritmo y exponencial. Esas funciones están implementado directamente en el hardware de hoy en día, y muy barato. O seleccionar juiciosamente algunos puntos y interpolar. O aproximado a través de splines. La magia rápida inversa de la raíz cuadrada del algoritmo (en realidad, su justificación) también pueden dar algunas ideas.

Sólo asegúrese de que esta operación es realmente relevante rendimiento sabia antes de ir por este camino.

Answer 3

Véase también stackoverflow optimizaciones-para-pow-con-const-no-entero-exponente de otros enfoques para el mismo problema.

Si nos fijamos en el gráfico de $x^{5/12}$ evidentemente tienes una forma diferente cerca de $x=0$, de cerca de $x=1$. Este es un buen indicio de que la simple polinomios no se ajustan tan bien a más de un dominio, incluyendo tanto $0$$1$.

Por lo que yo consideraría una reducción de la fortaleza:

considerar que en representación de coma flotante $x = m \cdot 2^k$ donde $0.5 \le m < 1$ así $$x^{5/12} = m^{5/12} \cdot 2^{5k/12}$$

MiniMaxApproximation[x^(5/12), {x, {1/2, 1}, 2, 2}]

Offers a rational O(2,2) polynomnial with a relative error less than $5.7\times10^{-7}$ for $m$ on $\left[\frac12,1\right]$.

For the second part of the strength reduction you can either use $\left[2,\frac {5k}{12}\right]$ or, since $k$ is an integer, optimise it with a 12 element lookup table and some modulo arithmetic.

Noting that, on most machines, division is very slow (and not good for numerical accuracy) I also considered simple polynomials and

MiniMaxApproximation[x^(5/12), {x, {1/2, 1}, 7, 0}]

Offers an O(7) polynomial with a relative error less than $2.4\times10^{-8}$ o mejor que la mitad de un ULP en precisión simple. Que debe ser de aproximadamente la misma velocidad que el S(2,2). La integridad, aquí es en Forma de Horner

0.246873 + 
   x (F954 + 
     x (-2.71714 + 
        x (3.86663 + 
          x (-3.85026 + x (2.47365 + (-0.919353 + 0.15006 x) x)))))