Aproximacion de funciones con polinomios de Chebyshev

https://cs.famaf.unc.edu.ar/~hoffmann/mn/practico06.html

La aproximación de Chebyshev (y cómo ahorrar dinero, ganar amigos, e influir sobre las personas)

Texto original de Jason M. Sachs (2012).

Bueno… tal vez sea una exageración. No creo que pueda recomendar nada que le ayudaría a ganar amigos. No es mi punto fuerte.

Pero voy a tratar de convencerte de porqué tendrías que conocer la aproximación de Chebyshev, que es un método para encontrar como uno puede acercarse lo más posible a calcular el resultado de una función matemática, con un mínimo de esfuerzo conceptual y de poder computacional.

Vamos a explorar dos casos de uso:

Amy tiene un microcontrolador 8-bit de poca potencia y necesita calcular $\sqrt{x}$ para tener una estimación precisa del valor cuadrático medio (o RMS) de la potencia disipada en una resistencia de detección de corriente. Encuentra que la implementación estandard de sqrt( ) toma mucho tiempo de CPU, y no tiene mucho espacio libre en memoria, y se pregunta si no hay otra manera de calcular raices cuadradas que use menos tiempo y memoria.
Bill está diseñando un detector de dióxido de azufre ( $SO_2$ ), con un sensor que emite un voltaje dependiente de la concentración del gas $SO_2$ , y un microcontrolador que convierte ese voltaje en una medida de la concentración de $SO_2$ que puede ser mostrada en pantalla. No tiene muchos problemas con las limitaciones de tiempo del CPU, pero tiene que inventar de alguna forma una manera de convertir el voltaje del sensor de gas en una concentración, porque el sensor es bastante no-linear.

En el caso de Amy, tiene una función que parece difícil de calcular. En el caso de bill, tiene una función desconocida lo cual debe llegar a modelizar. Estos son dos casos de evaluación de función con dos prioridades muy distintas. Los dos están en sistemas embedidos, lo que me lleva a la Hipótesis Central del Cálculo Matemático en Sistemas Enbedidos (HCdCMeSE? le hace falta un mejor nombre, y no me gusta llamarlo la Hipótesis Sachs):

Los únicos sistemas embedidos que requieren cálculos precisos en coma flotante son las calculadoras de escritorio.

Dejémolos decantar un minuto. Tu computadora es una cosa maravillosa. Puede calcular pi hasta centenas de decimales en unos milisegundos, o modelizar los patrones del flujo del aire alrededor de las palas de un motor. Los orígenes de la computadora han sido regularmente centrados alrededor del cálculo matemático, con un enorme pique de desarrollo cerca del fin de la segunda guerra mundial y durante la primera época de los transistores sde estado sólido. Si uno necesitaba una computadora en los años 1940 o 1950, era para hacer algun cálculo matemático preciso. La computadoras eran caras de tener y operar; si uno necesitaba una estimación rápida, había reglas de cálculo. En esa época también se publicaban libros enteros llenos de tablas de funciones de ingeniería, entonces si se necesitaba calcular funciones de Bessel hasta 5 decimales, uno podía hacerlo buscando en una tabla.

La computadora de escritorio usa una estrategia de cálculo llamada “precisión completa”: dada una representación en bits, hay solo una respuesta correcta a sin(x), y es el patrón de bits que produce el error más pequeño a partir del cálculo matemático exacto. Alguien escribió y comprobó librerias matemáticas para hacerlo, y desde ahí, es solo una llamada a una función.

Los sistemas embedidos, por otro lado, tienden a contener procesadores baratos que hacen cosas con medidas producidas por conversores analógico-digitales. Uno puede hacer el cálculo hasta todos los 53 bits de la doble precisión IEEE-754, pero si sus lecturas CAD, o el acondicionamiento de señal que las preceden, son solo precisas a 0,1% (10 bits), entonces los errores matemáticos van a ser minimizados por los errores de acondicionamiento de señal.

Cuando tenés un sistema que puede tolerar errores matemáticos, y tiene recursos limitados por los costos, entonces de repente la prioridad del cálculo matemático cambia. Ahora tenemos que decidir qué necesitamos realmente cuando queremos calcular algo.

Tomemos primero el problema de Amy: necesitamos calcular una raiz cuadrada.

Este planteo es solo un pantallazo de lo que se necesita. Hay dos decisiones mayores que tenemos que hacer de entrada:

¿sobre qué rango de $x$ tiene que ser calculado $\sqrt{x}$ ?
¿con qué precisión necesitamos tener $\sqrt{x}$ ? (¿y cómo varia con $x$ ?)

Las librerías generales con precisión completa proveen respuestas fáciles a estas preguntas: cualquier entrada tiene que ser aceptada, y para cada entrada hay una salida correcta, incluso si la salida correcta es un error (por ejemplo, $\sqrt{x}$ para entradas negativas en sistemas que no usan números complejos).

En sistemas embedidos, los requisitos de precisión y de rango son decisiones a tomar en conciencia.

La precisión es algo fácil de entender: hay algoritmos que toman atajos para calcular una expresión, y cuanta más imprecisión es tolerable, más atajos podemos tomar.

El requisito del rango es un poco más sútil. Obviamente, tenemos que soportar el rango esperado de entradas de un algoritmo, y tal vez queremos dejar algun margen. ¿Pero porqué sería más caro calcular $\sqrt{x}$ sobre un rango más grande? Si le dijera que me preocupa más el rango de un cálculo que su precisión, ¿me creería?

Usemos un ejemplo concreto. Amy piensa que necesita calcular $\sqrt{u}$ sobre el rango $u \in [0,2;5,0]$ . Su primer instinto es usar tablas de búsqueda.

Cómo no evaluar funciones, parte 1: tablas de búsqueda

El primer intento va a ser usar tablas de búsqueda. La idea es simple: tenés un arreglo de valores precomputados para tu función, y los usa como una tabla para evaluar la función. Miremos una solución para $\sqrt{x}$ sobre el rango $[0,2; 5,0]$ usando una tabla de 32 valores (acá estamos usando Python, supongamos que Amy hace algo equivalente en C):

import math
def maketable0(func, a, b, n):
    table = func([(x+0.5)/n*(b-a)+a for x in range(0,n)])
    def tablefunc(x):
        if x < a or x > b:
            raise "out of range error"
        i = math.floor((x-a)/(b-a)*n)
        if i == n:
            i = n-1
        return table[i]
    return tablefunc
sqrtapproxe = maketable0(np.sqrt, 0.2, 5.0, 32

Usar una tabla de busqueda involucra los dos procesos siguientes, que son bastante simples pero valen la pena ser examinados:

al momento de inicialización, dividir el rango de entradas de la función en N intervalos y evaluar la función dentro de esos intervalos (usualmente los intervalos son de tamaño igual y la evaluación de la función se hace en el punto medio de cada intervalo).
al momento de ejecución:
- manejar situaciones de fuera de rango (o tirar un error o devolver un valor correspondiente a la entrada máxima/mínima)
- calcular el índice de cual intervalo es el de la función
- buscar en la tabla el valor en ese índice

Así queda el error de apoximación de “e”:

Observá el resultado que parece una escalera y el grafo de error un serrucho. Eso es típico para una solución simple con tabla. El error máximo es gruesamente proporcional a la pendiente de la función dentro de un intervalo (el error es más largo para valores pequeños de $x$ donde $\sqrt{x}$ es más inclinado), y proporcional al tamaño del intervalo.

Si queremos mejorar la situación, la próxima mejora que podemos usar es la interpolación linear:

en el momento de inicialización, dividir el rango de entradas de la función en N intervalos y evaluar la función en los N+1 puntos de límite (usualmente los intervalos son de tamaño igual).
en el momento de ejecución:
- manejar situaciones de fuera de rango (o tirar un error o devolver un valor correspondiente a la entrada máxima/mínima)
- calcular el índice $k$ de cual intervalo es el de la función
- buscar en la tabla el valor en ese índice y en el siguiente
- hacer una interpolación lineal: y = table[k]*(1-u) + table[k+1]*u, donde $u$ es el número entre 0 y 1 que representa la posición de la entrada en el intervalo.

import math
def maketable1(func, a, b, n):
    table = func([float(x)/n*(b-a)+a for x in range(0,n+1)])
    def tablefunc(x):
        if x < a or x > b:
            raise "out of range error"
        ii = (x-a)/(b-a)*n
        i = math.floor(ii)
        if i == n:
            i = n-1
        u = ii-i
        return table[i]*(1-u) + table[i+1]*u
    return tablefunc
sqrtapproxf = maketable1(np.sqrt, 0.2, 5.0, 32)

Ahí va el grafo de error de “f”:

Mucho mejor – observá que el grafo de error contiene secciones parabólicas. Eso es típico de una tabla de búsqueda con interpolación. El error máximo es gruesamente proporcional a la curvatura (segunda derivada) de la función dentro de un intervalo, y proporcional al cuadrado del tamaño del intervalo. Si queremos reducir el error por un factor de 10, solo necesitamos 3 veces más entradas en la tabla de búsqueda.

Este es uno de los pocos casos donde una tabla de interpolación puede ser más apropiada que una evaluación polinomial: el error máximo es alrededor del medio del error del polinomio de grado 5 que usamos en la aproximación “d”.

Para una comparación en el tiempo de ejecución: una búsqueda en tabla sin interpolación es equivalente a una función cuadrática (grado 2), y una tabla con interpolación es equivalente a una función de cuarto grado. Es solo una estimación, depende realmente de cosas como si hacemos cálculos en coma flotante o no, o si la tabla tiene un tamaño que es una potencia de 2 (debería serlo, así se puede usar operación sobre bits para acelerar el cálculo), y qué tan largo se demora buscar en una tabla con el procesador que está usando.

Hay solo unas cuantas situaciones, en mi opinión, donde una tabla de búsqueda conviene:

Si no hace falta mucha precisión, una tabla sin interpolación puede ser suficiente.
Si una función es muy “no-polinomial” entonces una tabla interpolada puede lograr una precisión superior en comparación con una evaluación polinomial en tiempo equivalente. (Parece ser el caso para $\sqrt{x}$ sobre el rango $[0,2;5,0]$ . Vamos a definir “no-polinomial” dentro de poco; por ahora solo pensalo como una función que tiene muchas olas o curvas torcidas.)
Si una función es muy difícil de aproximar con un polinomio sobre un rango entero pero requiere una precisión muy alta, una opción es partir el rango en 4 u 8 subintervalos y usar una tabla de búsqueda para conseguir los coeficientes de cada subintervalo.

Las tablas de búsqueda no convienen cuando se requiere un tamaño de tabla tan grande que uno se debe preocupar por el uso de la memoria.

También debés observar que el error en la aproximación “f” no está regularmente distribuido: el error máximo ocurre en la parte baja del rango de entrada y es mucho menor sobre el resto del rango. Podés usar trucos para mejorar la distribución del error, como usar intervalos de tamaño variable, pero no hay un proceso general para hacerlo de manera eficiente en el momento de la ejecución.

Ahora, Amy piensa en buscar su libro de análisis matemático y la Serie de Taylor:

Cómo no evaluar funciones, parte 2: Serie de Taylor

Esta es la serie de Taylor para $\sqrt{1+x}$ (ecuación proveida por Wikipedia):

Pero hay una pequeña nota que dice que esto solo anda para $x$ entre $-1$ y $1$ , significando que no funciona para calcular $\sqrt{2}$ o la raiz cuadrada de cualquier cosa arriba de 2. Para $\sqrt{5}$ , tiene que usar un truco, a saber que $\sqrt{5} = 2 * \sqrt{5/4}$ . Entonces ella decide de limitar el rango de entradas de la serie de Taylor para cubrir el rango de $\sqrt{1+x}$ para $x\in [-0,8;+0,25]$ , y si $u = 1 + x > 1,25$ entonces Amy hará el truco de dividir-por-cuatro:

def sqrt1(x):
# returns sqrt(1+x) for x between -0.8 and +0.25
   return 1 + x/2 - x*x/8 + x*x*x/16
def sqrtapprox(u):
# valid for x between 0.2 to 5.0
   if (u >= 1.25):
      return sqrt1(u/4-1)*2
   else:
      return sqrt1(u-1)

¿Esto es suficiente? Hmm, pensémoslo un poco. Queremos conocer el error entre una función y su(s) aproximacion(es). Grafiquemos ese error. Ahi va una función que plotea los resultados:

import numpy as np
import matplotlib.pyplot as plt
def showplots(f,approxlist,a,b):
    x = np.linspace(a,b,1000)
    plt.figure(1)
    plt.subplot(211)
    plt.plot(x,f(x))
    vfuncs = [np.vectorize(approx) for approx in approxlist]
    for vf in vfuncs:
        plt.plot(x,vf(x))
    plt.xlim(a,b)
    plt.ylabel('f(x) and approximations fa(x)')
    plt.subplot(212)
    for vf in vfuncs:
        plt.plot(x,f(x)-vf(x))
    plt.xlim(a,b)
    plt.ylabel('error = f(x)-fa(x)')
    plt.xlabel('x')
    plt.show()

Los resultados de showplots(np.sqrt,[sqrtapprox],0.2,5.0) nos muestran el error entre la función sqrt de precisión completa y nuestra aproximación:

Puaj. Un error en el peor caso de $0,038$ . Intentemos con más términos de la serie de Taylor. Esto requiere entender como funcionan los coeficientes. Un poco de algebra con papel y lapicera muestra que los coeficientes desde el término cuadrático son:

$a_2 =$	$-\frac{1}{4\times 2} = -\frac{1}{8}$
$a_3 =$	$+\frac{3\times 1}{6\times 4\times 2} = \frac{1}{16}$
$a_4 =$	$-\frac{5\times 3\times 1}{8\times 6\times 4\times 2} = \frac{5}{128}$
$a_5 =$	$+\frac{7\times 5\times 3\times 1}{10\times 8\times 6\times 4\times 2} = \frac{7}{256}$
$a_6 =$	$-\frac{9\times 7\times 5\times 3\times 1}{12\times 10\times 8\times 6\times 4\times 2} = -\frac{21}{1024}$

Y sigue así. Entonces probemos de vuelta:

def sqrt1a(x):
# returns sqrt(1+x) for x between -0.8 and +0.25
   return 1 + x/2 - x*x/8 + x*x*x/16 - 5.0/128*x*x*x*x
def sqrtapproxa(u):
# valid for x between 0.2 to 5.0
   if (u >= 1.25):
      return sqrt1a(u/4-1)*2
   else:
      return sqrt1a(u-1)
def sqrt1b(x):
# returns sqrt(1+x) for x between -0.8 and +0.25
   return 1 + x/2 - x*x/8 + x*x*x/16 - 5.0/128*x*x*x*x + 7.0/256*x*x*x*x*x
def sqrtapproxb(u):
# valid for x between 0.2 to 5.0
   if (u >= 1.25):
      return sqrt1b(u/4-1)*2
   else:
      return sqrt1b(u-1)

Estos son los resultados para sqrtapproxa (hasta grado 4):

Y los de sqrtapproxb (hasta grado 5):

Puaj. Hemos agregado dos términos más y solo bajamos el error hasta $0,015$ .

Esto es lo que tenés que saber sobre la serie de Taylor para evaluar funciones. (O mejor, leelo y lo que no entendés, creeme).

Para obtener los coeficientes de una serie de Taylor para una función $f$ en una región alrededor de un punto $x=x_0$ , evaluás $f(x_0)$ y todas las derivadas de $f(x)$ en $x=x_0$ . No hace falta mirar a la función en otros puntos; el comportamiento en $x=x_0$ , de alguna manera, nos permite saber lo que es $f(x)$ en otras partes. Para ciertas funciones (como $\sqrt{1+x}$ ) la región de convergencia es finita; para otras como $sin(x)$ o $\exp^x$ , la región de convergencia es en todos lados.

La serie de Taylor funciona sobre la mayoría de las funciones cotidianas, porque esas funciones tienen una propiedad llamada “analicidad”: todas sus derivadas existen y son entonces continuas y lisas, y gracias a la magia de las matemáticas eso significa que mirar a una función analítica solo en 1 punto te da suficiente información para calcularla en alguna zona cercana.

Cerca de $x=x_0$ , la serie de Taylor es muy buena. El error entre la función exacta y los primeros $N$ términos de la serie de Taylor van a ser extremadamente pequeño alrededor de $x=x_0$ . A medida que te alejás de $x=x_0$ , el error empieza a crecer. En general, cuanto más términos usás, mejor será la aproximación alrededor de $x=x_0$ , pero luego la aproximación empieza a andar mal en el límite del rango, diverge más rápido.

Para resumir: la aproximación por serie de Taylor no distribuye el error de aproximación equitativamente. El error es pequeño en $x=x_0$ , pero ese error tiende a ser mucho más grande en los bordes de su rango.

Ahora le voy a mostrar dos intentos más para aproximar $\sqrt{x}$ para $x \in [0,2;5,0]$ , pero va a tener que esperar hasta que le explique como los conseguí.

def sqrt1c(x):
# valid for x between 0.2 to 1.25
   coeffs = [0.20569678, -0.94612133, 1.79170646, -1.89014568, 1.66083189, 0.17814197]
   y = 0
   for c in coeffs:
      y = y * x + c
   return y
def sqrtapproxc(u):
# valid for u between 0.2 to 5.0
   if (u >= 1.25):
      return sqrt1c(u/4)*2
   else:
      return sqrt1c(u)
def sqrtapproxd(x):
# valid for x between 0.2 to 5.0
   coeffs = [0.00117581,-0.01915684,0.12329254,-0.41444219,1.04368339, 0.26700714]
   y = 0
   for c in coeffs:
      y = y * x + c
   return y

La aproximación “c” parte el rango entre $[0,2; 1,25])$ y $[1.25; 5]$ como antes; la aproximación “d” opera sobre el rango entero $[0,2; 5]$ . Las dos son polinomios de grado 5. Este es el desempeño de la aproximación “c”:

Podés ver un pico de error de $0,00042$ . Eso es alrededor de 35 veces menos que el error de la aproximación “b” (la serie de Taylor con polinomios de grado 5 y la misma partición de rango).

Este el el desempeño de la aproximación “d”:

Hay un error en peor caso de $0,011$ , lo que es cerca de $30%$ menor al error de la aproximación “b”. ¡Mejoramos el error de peor caso, sin tener que partir el rango, con un polinomio de mismo grado!

Algunas cosas para observar en este ejercicio:

Las buenas funciones de aproximación tienen grafos de error con una “oleada” característica. Idealmente forman una función llamada “minimax”, donde los picos positivos y negativos se alinean todos. Eso reparte el error sobre todo el rango del argumento de entrada, y puede tener una precisión considerablemente mejor que un polinomio de una serie de Taylor del mismo grado. En práctica, es difícil conseguir el grafo de error minimax ideal, pero hay un método fácil que se acerca de ello, que voy a explicar en un ratito.
La reducción de rango puede mejorar considerablemente la precisión de la evaluación de una función. Ambas aproximaciones “c” y “d” usan un polinomio de grado 5, pero la aproximación “d” usa todo el rango de entada (con un ratio de 25:1 entre el máximo y el mínimo), mientras la aproximación “c” solo usa una parte del rango de entrada (con un ratio de 6.25:1), y “c” es 25 veces más precisa.

Polinomios de Chebyshev

La clave para usar polinomios para evaluar funciones, es no pensar en polinomios como siendo compuestos de combinaciones lineales de $1$ , $x$ , $x^2$ , $x^3$ , etc., pero como una combinación de Polinomios de Chebyshev $T_n(x)$ . Como se puede leer más en detalle en la página Wikipedia, estos tienen las propiedades siguientes:

son funciones minimax sobre el rango $[-1;+1]$ , es decir, todos sus mínimos/máximos son ±1 (ver imagen abajo)
son ortogonales sobre el rango $[-1;1]$ con peso $1/\sqrt{1-x^2}$
$T_0(x) = 1$ , $T_1(x) = x$ , $T_{n+1}(x) =2xT_n(x)-T_{n-1}(x)$
$T_n(x) = \cos(n \cos^{-1}(x))$

Supongamos que tenemos una función $f(x)$ sobre el rango $x\in[a;b]$ . Podemos expresar $f(x) = \sum c_k T_k(u)$ donde $u = \frac{2x-a-b}{b-a}$ y $x = \frac{b-a}{2}u + \frac{a+b}{2}$ , lo que es una transformación que mapea el intervalo $x\in[a;b]$ a $u\in[-1;1]$ . El problema entonces se transforma en encontrar los coeficientes $c_k$ , lo que puede hacerse usando los nodos de Chebyshev:

Son más fáciles de entender gráficamente; son simplemente proyecciones sobre el eje $x$ de puntos de un semicírculo dispuestos a igual distancia:

Podés ver que los nodos son regularmente distanciados cerca de $0$ y más comprimidos cerca de ± 1.

Para aproximar una función con una combinación linear de los primeros $N$ polinomios de Chebyshev ( $k=0$ a $N-1$ ), el coeficiente $c_k$ es simplemente igual a $A(k)$ multiplicado por el promedio de los productos $T_k(u)f(x)$ evaluados en los N nodos de Chebyshev, donde $A(0)=1$ y $A(k)=2$ para $k>1$ .

Ilustremos el proceso más claramente con un ejemplo.

Supongamos que $f(x) = \frac{1}{3} x^3 + 2x^2 + x - 10$ sobre el rango $[-1;3]$ , con $u = \frac{x-1}{2}, x = 2u+1$ :

Usemos $N=5$ , los nodos de Chebyshev para $N=5$ son $u=-0.951057, -0.587785, 0, +0.587785, +0.951057$ . Se corresponden con $x=-0.902113, -0.175571, 1, 2.175771, 2.902113$ , y la función $f(x)$ evaluada en esos nodos vale $y=-9.51921, -10.11572, -6.66667, 5.07419, 17.89408$ .

Para $c_0$ , dado que $T_0(u)f(x) = f(x)$ , calculamos el valor promedio de y = $-9.51921-10.11572-6.66667+5.07419+17.89408)/5 = -0.66667$

Para $c_1$ , dado que $T_1(u)f(x) = u\times f(x)$ , calculamos $2 \times$ el valor promedio de $u\times y =$ $2\times(-9.51921\times-0.951057+-10.11572\times-0.587785 +-6.66667\times 0+5.07419\times 0.587785+17.89408\times 0.951057) = 14$

Para $c_2$ , dado que $T_2(u)=2u T_1(u) - T_{n-1}(u)$ , calculamos $2 \times$ el valor promedio de $(2u^2 - 1) \times y \rightarrow c_2 = 6$

Para $c_3$ , calculamos $2 \times$ el valor promedio de $T_3(u) \times y = (4u^3 - 3u) \times y \rightarrow c_3 = 0.66667$

Para $c_4$ , calculamos $2 \times$ el valor promedio de $T_4(u) \times y = (8u^4 - 8u^2 + 1) \times y \rightarrow c_4 = 0$ .

La conclusión acá es que $f(x) = -\frac{2}{3}T_0(u) + 14T_1(u) + 6T_2(u) + \frac{2}{3}T_3(u)$ ; si procesás esto vamos a ver que es igual a $f(x) = \frac{1}{3}x^3 + 2x^2 + x - 10$ .

Entonces ¿cuál es el objetivo acá…? Bueno, si ya tenías los coeficientes de un polinomio, no hay punto en usar los polinomios de Chebyshev. (Si el rango de entrada es muy distinto de $[-1;1]$ , una transformación lineal de $x$ a $u$ donde $u\in [-1;1]$ te dará un polinomio $u$ que tiene mejor estabilidad numérica, y los polinomios de Chebyshev son una manera de hacer esa transformación. Si querés aproximar el polinomio con uno de menor grado, los polinomios de Chebyshev te darán la mejor aproximación.

Por ejemplo, usemos una función cuadrática para aproximar $f(x)$ sobre el rango de entrada. Todo lo que necesitamos hacer es truncar los coeficientes de Chebyshev. Calculemos:

$f_2(x) = -\frac{2}{3}T_0(u) + 14T_1(u) + 6T_2(u) = -\frac{2}{3} + 14u + 6\times(2u^2-1) = -\frac{20}{3} + 14u +12u^2 = -\frac{20}{3} + 7(x-1) + 3(x-1)^2 = 3x^2 + x - \frac{32}{3}.$

¡Listo! Tenemos una función cuadrática que es muy cercana a la función cúbica original, y la magnitud del error es solo el coeficiente del polinomio de Chebyshev que hemos sacado: $\frac{2}{3}$ .

Podés ver que los dos polinomios ( $\frac{1}{3}x^3 + 2x^2 + x - 10$ y $3x^2 + x - \frac{32}{3}$ ), expresados en su forma normal, no tienen una relación obvia uno con el otro, mientras que los coeficientes de Chebyshev son iguales salvo en el coeficiente $c_3$ que pusimos a cero para obtener una función cuadrática.

Desde mi punto de vista, eso es el uso principal de expresar una función en términos de sus coeficientes de Chebyshev: el coeficiente $c_0$ nos dice el valor promedio de la función sobre su rango de entrada, $c_1$ nos dice la “linearitud” de la función, $c_2$ nos dice su “cuadratitud”, $c_3$ su “cubiquitud”, etc.

Ahí va un ejemplo más típico y menos trivial. Supongamos que queremos aproximar $log_2 x$ sobre el rango $[1;2]$ .

Una decomposición en coeficientes de Chebyshev hasta el grado 6 nos da lo siguiente:

$c_0 =$	$+0.54311$
$c_1 =$	$+0.49505$
$c_2 =$	$-0.042469$
$c_3 =$	$+0.0048577$
$c_4 =$	$-6.2508\times 10^{-4}$
$c_5 =$	$+8.5757\times 10^{-5}$
$c_6 =$	$-1.1996\times 10^{-5}$

Podés ver que cada coeficiente es solo una fracción de magnitud del coeficiente anterior. Eso es un indicador de que la función en cuestión (sobre el rango de entrada especificado) se presta bien a una aproximación polinomial. Si los coeficientes de Chebyshev no bajan fuertemente alrededor del coeficiente 5, yo diría que esa función es “no-polinomial” sobre el rango de interés.

El error máximo para la aproximación de Chebyshev de grado 6 de $log_2 x$ sobre $x\in[1;2]$ es tan solo $2,2 \times 10^{-6}$ , lo que es aproximadamente el valor de coeficiente de Chebyshev siguiente.

Si truncamos la serie a una aproximación de Chebyshev de grado 4, obtenemos un error máximo que es casi igual a $c_5 = 8,6\times 10^{-5}$ :

¿No está mal, si?

Más ejemplos de funciones elementales

Ahí va una lista de funciones comunes, con sus primeros 6 coeficientes de Chebyshev:

$f(x)$	range	c₀	c₁	c₂	c₃	c₄	c₅
$4\sin \pi x$	[-0.5,0.5]	0	1.1336	0	-0.13807	0	0.0045584
$\sin \pi x$	[-0.25,0.25]	0	0.72638	0	-0.01942	0	0.00015225
$\cos \pi x$	[-0.5,0.5]	0.472	0	-0.4994	0	0.027985	0
$\cos \pi x$	[-0.25,0.25]	0.85163	0	-0.14644	0	0.0019214	0
$\sqrt{x}$	[1,4]	1.542	0.49296	-0.040488	0.0066968	-0.0013836	0.00030211
$\log_2 x$	[1,2]	0.54311	0.49505	-0.042469	0.0048576	-0.00062481	8.3994e-05
$e^x$	[0,1]	1.7534	0.85039	0.10521	0.0087221	0.00054344	2.7075e-05
$\frac{2}{\pi}\tan^{-1} x$	[-1,1]	0	0.5274	0	-0.030213	0	0.0034855
$\frac{1}{1 + e^{-x}}$	[-1,1]	0.5	0.23557	0	-0.0046202	0	0.00011249
$\frac{1}{1 + e^{-x}}$	[-3,3]	0.5	0.50547	0	-0.061348	0	0.01109
$\frac{1}{1 + x^2}$	[-1,1]	0.70707	0	-0.24242	0	0.040404	0
$\frac{1}{1 + x^2}$	[-3,3]	0.30404	0	-0.29876	0	0.12222	0

Otros comentarios

Una buena referencia sobre las funciones de Chebyshev es Numerical Recipes de Press, Teukolsky, Vetterlingm y Flannery, que cubre la aproximación de Chebyshev en detalle.

Hay unas cosas para saber cuando evaluamos funciones de Chebyshev:

es mejor evaluar las funciones directamente que tratar de convertir las aproximaciones de Chebyshev en una forma polinomial estandard (es decir, si querés calcular $3T_0(x) + 5T_1(x) - 0.1T_2(x)$ , no lo convirtas a un polinomio pero calculá los valores de $T$ directamente.) Numerical Recipes recomienda la formula de recurrencia de Clenshaw, pero yo uso simplemente el algoritmo siguiente:
```
u = (2*x-a-b)/(b-a)
Tprev = 1
T = u
y = c[0]
for i = 1:n
   y=y+T*c[i]
   Tnext = 2*u*T - Tprev
   Tprev = T
   T = Tnext
```
recordá de convertir $x$ en $u$ , eso mantiene la estabilidad numérica del algoritmo, porque los polinomios tienden a ser bien acondicionados para evaluar sus argumentos entre -1 y +1, y mal acondicionados con valores muy grandes o muy pequeños.

Ahí va una clase en Python para calcular y aplicar coeficientes de Chebyshev:

import math
import numpy as np
def chebspace(npts):
    t = (np.array(range(0,npts)) + 0.5) / npts
    return -np.cos(t*math.pi)
def chebmat(u, N):
    T = np.column_stack((np.ones(len(u)), u))
    for n in range(2,N+1):
        Tnext = 2*u*T[:,n-1] - T[:,n-2]
        T = np.column_stack((T,Tnext))
    return T
class Cheby(object):
    def __init__(self, a, b, *coeffs):
        self.c = (a+b)/2.0
        self.m = (b-a)/2.0
        self.coeffs = np.array(coeffs, ndmin=1)
    def rangestart(self):
        return self.c-self.m
    def rangeend(self):
        return self.c+self.m
    def range(self):
        return (self.rangestart(), self.rangeend())
    def degree(self):
        return len(self.coeffs)-1
    def truncate(self, n):
        return Cheby(self.rangestart(), self.rangeend(), *self.coeffs[0:n+1])
    def asTaylor(self, x0=0, m0=1.0):
        n = self.degree()+1
        Tprev = np.zeros(n)
        T     = np.zeros(n)
        Tprev[0] = 1
        T[1]  = 1
        # evaluate y = Chebyshev functions as polynomials in u
        y     = self.coeffs[0] * Tprev
        for co in self.coeffs[1:]:
            y = y + T*co
            xT = np.roll(T, 1)
            xT[0] = 0
            Tnext = 2*xT - Tprev
            Tprev = T
            T = Tnext
        # now evaluate y2 = polynomials in x
        P     = np.zeros(n)
        y2    = np.zeros(n)
        P[0]  = 1
        k0 = -self.c/self.m
        k1 = 1.0/self.m
        k0 = k0 + k1*x0
        k1 = k1/m0
        for yi in y:
            y2 = y2 + P*yi
            Pnext = np.roll(P, 1)*k1
            Pnext[0] = 0
            P = Pnext + k0*P
        return y2
    def __call__(self, x):
        xa = np.array(x, copy=False, ndmin=1)
        u = np.array((xa-self.c)/self.m)
        Tprev = np.ones(len(u))
        y = self.coeffs[0] * Tprev
        if self.degree() > 0:
            y = y + u*self.coeffs[1]
            T = u
        for n in range(2,self.degree()+1):
            Tnext = 2*u*T - Tprev
            Tprev = T
            T = Tnext
            y = y + T*self.coeffs[n]
        return y
    def __repr__(self):
        return "Cheby%s" % (self.range()+tuple(c for c in self.coeffs)).__repr__()
    @staticmethod
    def fit(func, a, b, degree):
        N = degree+1
        u = chebspace(N)
        x = (u*(b-a) + (b+a))/2.0
        y = func(x)
        T = chebmat(u, N=degree)
        c = 2.0/N * np.dot(y,T)
        c[0] = c[0]/2
        return Cheby(a,b,*c)

Podés guardar esto en un archivo cheby.py y luego usarlo como sigue (acá hacemos caber $sin x$ entre 0 y $\frac{\pi}{2}$ con un polinomio de quinto grado, y calculamos nuestra aproximación en los ángulos $0, \frac{\pi}{6}, \frac{\pi}{4} y \frac{\pi}{3}$ ):

>>> import cheby
>>> import numpy as np
>>> import math
>>> c = cheby.Cheby.fit(np.sin,0,math.pi/2,5)
>>> c
Cheby(0.0, 1.5707963267948966, 0.60219470125550711, 0.51362516668030367, -0.10354634422944738, -0.013732035086651754, 0.001358650338492214, 0.00010765948465629727)
>>> c(math.pi*np.array([0, 1.0/6, 1.0/4, 1.0/3]))
array([  6.21628624e-06,   5.00003074e-01,   7.07099696e-01,
         8.66028717e-01])

Aproximación de funciones para datos empíricos

Vuelvamos al principio de nuestra charla, al problema de Bill. Tiene un sensor no-lineal relacionando una cantidad física (concentración de gas) con una medición de voltaje. El problema que tiene es como invertir esa relación para que pueda estimar la cantidad física original: $y = f(x)$ donde $x$ es la medición e $y$ es la cantidad física.

La aproximación de Chebyshev funciona bien en este caso también, pero tenemos un problema. No podemos simplemente evaluar $f(x)$ en los nodos de Chebyshev para obtener los coeficientes de Chebyshev, porque no estamos seguros de cuánto vale $f(x)$ en realidad. En lugar de eso, tenemops que hacer unas mediciones de referencia para determinar esos coeficientes, con cantidades físicas conocidas (es decir, concentraciones de gas), y hacer una tabla con esos valores x e y, y luego usarlos para estimar los coeficientes. Tenemos dos elecciones posibles para seleccionar esas mediciones:

tratar de planificar esas mediciones para que los valores x estén cerca de los nodos de Chebyshev sobre el rango esperado. Eso nos permitiría usar menos mediciones.
Hacer muchas mediciones para cubrir bien todo el rango de entradas.

En ciertas situaciones, es más rápido y fácil hacer mediciones con referencias arbitrarias> (Por ejemplo, si está tratando de calibrar un micrófono, puede cambiar la amplitud y frecuencia de las señales sonoras de referencia automáticamente con equipamiento electrónico), y la segunda opción es viable. En otras situaciones, como con la concentración de gas, puede ser inconveniente o caro hacer muchas mediciones, y tendrá que seleccionar los valores de referencia con cuidado.

La mejor manera que encontré para estimar coeficientes de Chebyshev a partir de datos empíricos, es usando los mínimos cuadrados ponderados:

colectar las mediciones x y la entrada física y en vectores X e Y
para el vector de mediciones x, calcular u a partir de x, y calcular los valores de los polinomios de Chebyshev $T_0(u), T_1(u), \ldots, T_n(u)$ y guardar cada uno de esos como vectores en una matriz T. Cada polinomio va a ser usado como una función base.
Calcular un valor de ponderación $w(x)$ para compensar la densidad de los valores $x$ (por ejemplo, si hay muchos valores x cerca de $x=0$ pero muy pocos cerca de $x=1$ , puede elegir un $w(x)$ bajo cerca de 0 y más alto cerca de !.) Expresar esos valores como una matriz diagonal W.
Calcular un vector C de coeficientes $c_0$ , $c_1$ , $c_2$ , $\ldots$ , $c_n$ usando los mínimos cuadrados ponderados para resolver C: $(T^TWT)C=T^TWY$ Para hacer eso, tendra que usar alguna librería de matrices. En Matlab se calcularía con C = (T'*diag(W)*T)\(T'*diag(W)*Y). En Python con Numpy se puede saltar la etapa de calcular la matriz T y directamente usar la función numpy.polynomial.chebyshev.chebfit sobre las mediciones x e y; el vector de ponderaciones es un parámetro opcional.

Yo no usaría un grado mayor a 5. Si encontrás que todavía no tenés una buena aproximación con un polinomio de grado 5, probablemente no está usando el método correcto; la función es demasiado puntiaguda u ondulada o cuspidosa. Los polinomios de grado mayor pueden tener problemas con la precisión numérica.

En todos casos, la aproximación funcional de datos empíricos es un poco más complicada que la aproximación de funciones conocidas. Asegúrese siempre de verificar dos veces su aproximación haciendo un grafo de los datos originales y de la función que encontró.

Resumen

La proxima vez que tenés que usar la aproximación de funciones, ¡probá la aproximación de Chebyshev! No solo es probablemente la mejor y más fácil manera de aproximar una función con un polinomio, pero también te va a indicar qué tan bien los polinomios aproximan la función en cuestión, por el comportamiento de los coeficientes de Chebyshev.

Asegurate de elegir la precisión y el rango de entradas de manera razonable. Una precisión demasiado grande y un rango demasiado ancho va a requerir polinomios de grado más alto y tiempos de ejecución mayores, y los polinomios de alto grado son a menudo mal acondicionados numéricamente.

Una evaluación de funciones más eficiente significa que podrás usar menos recursos del CPU para una precisión deseada. Entonces esto puede hacerte ahorrar dinero e influir sobre tus colegas.

¡Feliz aproximación!

Diseño Electrónico

sábado, 19 de septiembre de 2020