¿Cómo puedo desarrollar rutinas numéricas para la evaluación de mis propias funciones especiales?

Question

Esta pregunta ha sido enviada a ComputationalScience.SE aquí .

Cuando realizo un trabajo computacional, a menudo me encuentro con una función univariante, definida en términos de una ecuación integral o diferencial, que me gustaría evaluar rápidamente (digamos, millones de veces por segundo) en un intervalo específico con una precisión determinada (digamos, una parte en $10^{10}$ ). Por ejemplo, la función $$ f(\alpha) = \int_{k=0}^\infty \frac{e^{-\alpha^2 k^2}}{k+1}\ \mathrm{d}k $$ en el intervalo $\alpha \in (0,10)$ surgió en un proyecto reciente. Ahora resulta que esta integral puede ser evaluada en términos de funciones especiales estándar (en particular, $\operatorname{Ei}(z)$ y $\operatorname{erfi}(z)$ ), pero supongamos que tenemos una función mucho más complicada para la que no se conoce dicha evaluación. ¿Existe una técnica sistemática ¿Puedo aplicar para desarrollar mis propias rutinas numéricas para la evaluación de tales funciones?

Estoy seguro de que debe haber muchas técnicas por ahí, ya que parece que existen algoritmos rápidos para básicamente todas las funciones especiales comunes. Sin embargo, insisto en que el tipo de técnica que busco no debería depender de que la función tenga alguna estructura particular (por ejemplo, relaciones de recurrencia como $\Gamma(n+1) = n\Gamma(n)$ o fórmulas de reflexión como $\Gamma(z) \Gamma(1-z) = \pi \csc(\pi z)$ ). Idealmente, esta técnica funcionaría para casi cualquier función (lo suficientemente bien comportada) con la que me encuentre.

Puedes dar por sentado que tengo algo de lento método de evaluación de la función deseada (por ejemplo, integración numérica directa) con cualquier precisión, y que estoy dispuesto a hacer mucho trabajo de preprocesamiento con el método lento para desarrollar un método rápido.

Answer 1

Desgraciadamente, no existe un enfoque único que conduzca a implementaciones robustas, precisas y de alto rendimiento en todo el amplio universo de funciones especiales. A menudo hay que utilizar dos o más métodos para diferentes partes del dominio de entrada, y el trabajo de investigación e implementación necesario puede llevar semanas para las funciones elementales y meses para las funciones trascendentales superiores.

Dado que se requieren considerables conocimientos matemáticos y de programación para producir implementaciones de alta calidad, mi primera recomendación sería utilizar y aprovechar las bibliotecas matemáticas existentes en la medida de lo posible. Pueden ser bibliotecas comerciales, como la Bibliotecas numéricas NAG o el Bibliotecas numéricas IMSL de RogueWave o bibliotecas de código abierto como el Biblioteca científica de GNU (GSL) o el partes matemáticas y numéricas de la biblioteca Boost . También puede encontrar el código fuente pertinente en repositorios en línea como Algoritmos recogidos por Netlib en ACM TOMS .

En términos prácticos, en los procesadores modernos mejorados con SIMD, el uso extensivo de tablas ya no es aconsejable y las aproximaciones polinómicas (a trozos) suelen ser las más atractivas. La razón por la que los métodos basados en tablas han caído en desgracia en las arquitecturas de procesadores de alto rendimiento es que en la última década el rendimiento de las unidades funcionales (medido en FLOPS) ha crecido mucho más rápido que el rendimiento de los subsistemas de memoria (medido en GB/seg). El razonamiento del siguiente documento coincide con mi propia experiencia profesional:

Marat Dukhan y Richard Vuduc, "Methods for high-throughput computation of elementary functions". En Procesamiento paralelo y matemáticas aplicadas , pp. 86-95. Springer, 2014. ( diapositivas )

En términos de rendimiento, las aproximaciones polinómicas se benefician de la operación de multiplicación- suma fusionada (FMA) presente en el hardware de los procesadores modernos (tanto CPU como GPU). Esta operación también ayuda a reducir los errores de redondeo y ofrece cierta protección contra la cancelación sustractiva. Para obtener el menor error y la mejor eficiencia, se recomienda utilizar aproximación minimax .

Herramientas de uso común como Arce y Mathematica cuentan con facilidades incorporadas para generarlas. Aunque generan aproximaciones que son (muy cercanas a) óptimas en el sentido matemático, no suelen tener en cuenta el error en el que se incurre por la representación de los coeficientes y la evaluación de las operaciones en una precisión de punto flotante limitada. La página web Sollya ofrece esta funcionalidad a través de su fpminimax de la empresa. Por último, puede escribir su propio código de aproximación, que probablemente se basará en el Algoritmo Remez .

Para algunas funciones, las aproximaciones polinómicas no son realmente prácticas, ya que se necesitarían demasiados términos para alcanzar la doble precisión IEEE-754. En esos casos, se puede optar por una de las dos estrategias.

La primera estrategia consiste en transformar inteligentemente los argumentos de entrada, utilizando aritmética básica y funciones elementales sencillas, de manera que la función resultante se "comporte bien" con respecto a la aproximación polinómica. Por lo general, dicha transformación tiende a "linealizar" la función a aproximar. Un buen ejemplo didáctico de este enfoque es el cálculo de erfc en el siguiente documento:

M. M. Shepherd y J. G. Laframboise, "Chebyshev Approximation of $(1 + 2x)\exp(x^2)\operatorname{erfc} x$ en $0 \leqslant x < \infty$ ". Matemáticas de la computación Vol. 36, nº 153 (enero de 1981), pp. 249-253 (en línea)

El segundo enfoque consiste en utilizar el cociente de dos polinomios, es decir, una aproximación racional, por ejemplo en forma de Aproximador de Padé . Las herramientas ya mencionadas pueden ayudar a ello; también hay una copiosa cantidad de literatura publicada sobre la aproximación racional, que en general es un problema más difícil que la aproximación polinómica.

En el caso de las funciones especiales (a diferencia de las funciones elementales), las aproximaciones polinómicas y racionales directas suelen ser inexactas y/o ineficaces. Requieren la aplicación de conceptos matemáticos más avanzados, como expansiones asintóticas, relaciones de recurrencia y expansiones de fracciones continuas. Incluso si el uso de estos conceptos resuelve el problema matemáticamente, puede haber problemas numéricos, por ejemplo cuando se evalúan fracciones continuas en la dirección de avance. No es de extrañar que se hayan escrito libros enteros sobre la evaluación por ordenador de ciertas funciones, como las funciones de Bessel y las funciones de Mathieu.

A continuación, doy una rápida visión de la literatura útil, empezando por la cobertura de los fundamentos matemáticos, pasando por los métodos adecuados para las funciones elementales y las funciones especiales simples como erfc y tgamma y, por último, métodos avanzados para funciones especiales que son más difíciles de calcular tanto en términos de rendimiento como de precisión. Obviamente, esto sólo puede arañar la superficie, ya que se puede encontrar mucho material relevante sobre funciones particulares en artículos individuales, por ejemplo, en revistas y actas de AMS, SIAM, ACM e IEEE.

Gran parte de la bibliografía aún no se ha puesto al día con los modernos entornos de hardware y software, en particular con la presencia de la operación FMA y las arquitecturas SIMD. En cuanto a los códigos informáticos robustos para la evaluación de funciones matemáticas, cabría desear una cooperación más estrecha entre las matemáticas y la ciencia, por un lado, y la informática y la ingeniería informática, por otro. Entre los trabajos que se presentan a continuación, los de Markstein y Muller son los más avanzados en este sentido.

Milton Abramowitz e Irene A. Stegun (eds.), "Handbook of Mathematical Functions. With Formulas, Graphs, and Mathematical Tables". Nueva York, NY: Dover 1972 ( versión en línea )

Frank Olver, et. al. (eds.), "NIST Handbook of Mathematical Functions". Nueva York, NY: Cambridge University Press 2010 ( versión en línea )

A. Erdelyi, et. al., "Higher Transcendental Functions". Vol. 1-3. Nueva York, NY: McGraw-Hill 1955

Oskar Perron, "Die Lehre von den Kettenbrüchen, 3ª edición". Vol. 1+2. Stuttgart (Alemania): Teubner 1954, 1957

---

John F. Hart, "Computer Approximations". Malabar, FL: Krieger Publishing 1978

William J. Cody y William Waite, "Software Manual for the Elementary Functions". Englewood Cliffs, NJ: Prentice-Hall 1980

Peter Markstein, "IA-64 y funciones elementales". Upper Saddle River, NJ: Prentice-Hall 2000

Jean-Michel Muller, "Funciones elementales. Algorithms and Implementation 3rd. ed.". Birkhäuser 2016

Nelson H. F. Beebe, "The Mathematical-Function Computation Handbook". Springer 2017

Jean-Michel Muller, et. al., "Handbook of Floating-Point Arithmetic 2nd ed.". Birkhäuser 2018

---

Nico M. Temme, "Funciones especiales. Una introducción a las funciones clásicas de la física matemática". Nueva York, NY: Wiley 1996

Amparo Gil, Javier Segura, y Nico M. Temme, "Numerical Methods for Special Functions". SIAM 2007

---

Frank W. J. Olver, "Asymptotics and Special Functions". Natick, MA: A K Peters 1997

Jet Wimp, "Computation with Recurrence Relations". Boston, MA: Pitman 1984

---

A. N. Khovanskii, "The application of continued fractions and their generalizations to problems in approximation theory". Groningen (Países Bajos): Noordhoff 1963

A. Cuyt, et. al., "Handbook of Continued Fractions for Special Functions". Springer 2008