Para datos ruidosos o de estructura fina, ¿hay mejores cuadraturas que la regla del punto medio?

Solo las dos primeras secciones de esta larga pregunta son esenciales. Los otros son solo para ilustración.

Antecedentes

Las cuadraturas avanzadas como Newton-Cotes, Gauß-Legendre y Romberg compuestos de mayor grado parecen estar destinadas principalmente a casos en los que se puede muestrear finamente la función pero no integrarse analíticamente. Sin embargo, para funciones con estructuras más finas que el intervalo de muestreo (ver Apéndice A para un ejemplo) o medición de ruido, no pueden competir con enfoques simples como el punto medio o la regla trapezoidal (ver Apéndice B para una demostración).

Esto es algo intuitivo ya que, por ejemplo, la regla compuesta de Simpson esencialmente "descarta" una cuarta parte de la información al asignarle un peso menor. La única razón por la cual tales cuadraturas son mejores para funciones suficientemente aburridas es que el manejo adecuado de los efectos de borde supera el efecto de la información descartada. Desde otro punto de vista, para mí es intuitivamente claro que para funciones con una estructura fina o ruido, las muestras que están alejadas de los límites del dominio de integración deben ser casi equidistantes y tener casi el mismo peso (para un gran número de muestras ) Por otro lado, la cuadratura de tales funciones puede beneficiarse de un mejor manejo de los efectos de borde (que para el método del punto medio).

Pregunta

Suponga que deseo integrar numéricamente datos unidimensionales ruidosos o bien estructurados.

El número de puntos de muestreo es fijo (debido a que la evaluación de la función es costosa), pero puedo ubicarlos libremente. Sin embargo, yo (o el método) no puedo colocar puntos de muestreo de forma interactiva, es decir, en función de los resultados de otros puntos de muestreo. Tampoco sé de antemano posibles regiones problemáticas. Entonces, algo como Gauß – Legendre (puntos de muestreo no equidistantes) está bien; la cuadratura adaptativa no lo es, ya que requiere puntos de muestreo colocados interactivamente.

• ¿Se ha sugerido algún método que vaya más allá del método del punto medio para tal caso?

• O: ¿Hay alguna prueba de que el método del punto medio es mejor en tales condiciones?

• Más en general: ¿Hay algún trabajo existente sobre este problema?

Apéndice A: ejemplo específico de una función bien estructurada

Deseo estimar para: cony. Una función típica se ve así:$\int_0^1f(t)\, \mathrm{d}t$

$$f(t) = \sum_{i=1}^{k} \frac{\sin(ω_i t-φ_i)}{ω_i},$$ $φ_i∈ [0,2π]$$\log{ω_i} ∈ [1,1000]$

Elegí esta función para las siguientes propiedades:

• Se puede integrar analíticamente para un resultado de control.
• Tiene una estructura fina en un nivel que hace que sea imposible capturarlo todo con la cantidad de muestras que estoy usando ( ).$<10^2$
• No está dominado por su fina estructura.

Apéndice B: Benchmark

Para completar, aquí hay un punto de referencia en Python:

import numpy as np
from numpy.random import uniform
from scipy.integrate import simps, trapz, romb, fixed_quad

begin = 0
end   = 1

def generate_f(k,low_freq,high_freq):
    ω = 2**uniform(np.log2(low_freq),np.log2(high_freq),k)
    φ = uniform(0,2*np.pi,k)
    g = lambda t,ω,φ: np.sin(ω*t-φ)/ω
    G = lambda t,ω,φ: np.cos(ω*t-φ)/ω**2
    f = lambda t: sum( g(t,ω[i],φ[i]) for i in range(k) )
    control = sum( G(begin,ω[i],φ[i])-G(end,ω[i],φ[i]) for i in range(k) )
    return control,f

def midpoint(f,n):
    midpoints = np.linspace(begin,end,2*n+1)[1::2]
    assert len(midpoints)==n
    return np.mean(f(midpoints))*(n-1)

def evaluate(n,control,f):
    """
    returns the relative errors when integrating f with n evaluations
    for several numerical integration methods.
    """
    times = np.linspace(begin,end,n)
    values = f(times)
    results = [
            midpoint(f,n),
            trapz(values),
            simps(values),
            romb (values),
            fixed_quad(f,begin,end,n=n)[0]*(n-1),
        ]

    return [
            abs((result/(n-1)-control)/control)
            for result in results
        ]

method_names = ["midpoint","trapezoid","Simpson","Romberg","Gauß–Legendre"]

def med(data):
    medians = np.median(np.vstack(data),axis=0)
    for median,name in zip(medians,method_names):
        print(f"{median:.3e}   {name}")

print("superimposed sines")
med(evaluate(33,*generate_f(10,1,1000)) for _ in range(100000))

print("superimposed low-frequency sines (control)")
med(evaluate(33,*generate_f(10,0.5,1.5)) for _ in range(100000))

(Aquí uso la mediana para reducir la influencia de los valores atípicos debido a funciones que solo tienen contenido de alta frecuencia. Por lo demás, los resultados son similares).

Las medianas de los errores de integración relativos son:

superimposed sines
6.301e-04   midpoint
8.984e-04   trapezoid
1.158e-03   Simpson
1.537e-03   Romberg
1.862e-03   Gauß–Legendre

superimposed low-frequency sines (control)
2.790e-05   midpoint
5.933e-05   trapezoid
5.107e-09   Simpson
3.573e-16   Romberg
3.659e-16   Gauß–Legendre

Nota: Después de dos meses y una recompensa sin resultado, publiqué esto en MathOverflow .

En primer lugar, creo que no entiendes el concepto de cuadratura adaptativa. La cuadratura adaptativa no implica "colocar puntos de muestra de forma interactiva". La idea general detrás de la cuadratura adaptativa es idear un esquema que integrará una determinada función a un cierto error absoluto o relativo (estimado) con la menor cantidad posible de evaluaciones de funciones.

Una segunda observación: escribe "El número de puntos de muestreo es fijo (debido a que la evaluación de la función es costosa), pero puedo ubicarlos libremente". Creo que la idea debería ser que el número de puntos de muestreo (o evaluaciones de funciones en terminología en cuadratura) debería ser lo más pequeño posible (es decir, no fijo).

Entonces, ¿cuál es la idea detrás de la cuadratura adaptativa implementada en QUADPACK, por ejemplo?

1. El ingrediente básico es una regla de cuadratura "anidada": esta es una combinación de dos reglas de cuadratura donde una tiene un orden (o precisión) más alto que la otra. ¿Por qué? Basado en la diferencia entre estas reglas, el algoritmo puede estimar el error de cuadratura (por supuesto, el algoritmo usará el más preciso como resultado de referencia). Ejemplos podrían ser la regla trapezoidal con nodos y nodos. En el caso de QUADPACK, las reglas son reglas de Gauss-Kronrod. Estas son reglas de cuadratura interpolatoria que usan una regla de cuadratura de Gauss-Legendre de cierto orden$2^{n}$$2^{n+1}$$N$y una extensión óptima de esta regla. Esto significa que se puede obtener un orden de cuadratura más alto reutilizando los nodos de Gauss-Legendre (es decir, las evaluaciones de funciones costosas) con diferentes pesos y agregando varios nodos adicionales. En otras palabras, la regla original de orden de Gauss-Legendre integrará todos los polinomios de grado exactamente mientras que la regla extendida de Gauss-Kronrod integrará algún polinomio de orden superior exactamente. Una regla clásica es el G7K15 (Gauss-Legendre de 7 ° orden con Gauss-Kronrod de 15 ° orden). La magia es que los 7 nodos de Gauss-Legendre son un subconjunto de los 15 nodos de Gauss-Kronrod, así que con 15 evaluaciones de funciones, ¡tengo una evaluación en cuadratura junto con una estimación de error!$N$$2N-1$

2. El siguiente ingrediente es una estrategia de "divide y vencerás". Suponga que suelta este G7K15 en su integrando y observa un error de cuadratura que, según su gusto, es demasiado grande. QUADPACK subdividirá el intervalo original en dos subintervalos igualmente espaciados. Y luego volverá a evaluar las dos subintegrales usando la regla básica, G7K15. Ahora, el algoritmo tiene una estimación de error global (que debería ser más baja que la primera) pero también dos estimaciones de error locales. Elige el intervalo con el error más grande y divide este en dos. Se estiman dos nuevas integrales y se actualiza el error global. Y así sucesivamente hasta que el error global esté por debajo de su objetivo solicitado o se haya superado el número máximo de subdivisiones.

Por lo tanto, le desafío a actualizar su código anterior utilizando el scipy.quadmétodo Quizás, en el caso de un integrando con mucha "estructura fina", necesite aumentar el número máximo de subdivisiones (la limitopción). También puedes jugar con los parámetros epsabsy / o epsrel.

Sin embargo, si solo tiene datos experimentales, veo dos posibilidades.

1. Si tiene la oportunidad de seleccionar los puntos de medición, es decir, los valores de , los seleccionaría equidistantemente y preferiblemente como una potencia de para que pueda aplicar una regla trapezoidal anidada (y aprovechar la extrapolación de Romberg).$t$$2$
2. Si no tiene medios para elegir los nodos, es decir, las mediciones se realizan en momentos aleatorios, la mejor opción en mi opinión sigue siendo la regla del trapecio.

No estoy convencido de que su código demuestre algo fundamental sobre las diversas reglas de cuadratura y qué tan bien lo hacen contra el ruido y la estructura fina, y creo que es probable que si elige varias estructuras de multas diferentes encuentre algo diferente. Aquí está el teorema:

Ningún método de cuadratura puede dar un error absoluto o relativo bajo contra una función con variación total ilimitada. En un sistema de coma flotante con redondeo unitario , tenemos la estimación$\mu$ donde es la suma en cuadratura que actúa sobre la implementación numérica de .

$$\left| \int_{a}^{b} f \, \mathrm{d}x - \hat{Q}[\hat{f}] \right| \le \left| \int_{a}^{b} f \, \mathrm{d}x - Q[f] \right| + \mu\left[ 4\int_{a}^{b} |f| \, \mathrm{d}x + \int_{a}^{b} |xf'| \, \mathrm{d}x \right]$$ $\hat{Q}$$\hat{f}$$f$

Prueba: deje que los nodos en cuadratura sean y los pesos en cuadratura (no negativos) sean y denotan sus aproximaciones de punto flotante por y . Suponga que satisface donde donde es el redondeo de la unidad. Luego $\{x_i\}_{i=0}^{n-1}$$\{w_i\}_{i=0}^{n-1}$$\hat{w}_{i}$$\hat{x}_i$$\hat{f}$$\hat{f}(x) = f(x)(1+2\delta)$$|\delta| \le \mu$$\mu$

$$\begin{align*} \hat{Q}[\hat{f}] &= \sum_{i=0}^{n-1} \hat{w}_i \otimes \hat{f}(\hat{x}_i) \\ &= \sum_{i=0}^{n-1} w_i (1+\delta^w_i)f(x_i + \delta_i^x x_i)(1+2\delta_i^{f})(1+\delta_i^{*}) \\ &\approx \sum_{i=0}^{n-1} w_i \left[f(x_i)+ \delta_i^x x_i f'(x_i) \right] (1+\delta^w_i + 2\delta_i^{f} + \delta_i^*) \\ &\approx \sum_{i=0}^{n-1} w_i f(x_i) + \sum_{i=0}^{n-1}\delta_i^x w_i x_i f'(x_i) + w_i f(x_i) (\delta^w_i + 2\delta_i^{f} + \delta_i^*) \\ \end{align*}$$ para que Esto supone que la suma se calcula sin error; multiplique por para descartar esa suposición. $$\begin{align*} |\hat{Q}[\hat{f}] - Q[f]| &\le \mu \sum_{i=0}^{n-1}w_i(|x_i f'(x_i)| + 4|f(x_i)|) \\ &\approx 4\mu \int |f| \, \mathrm{d}x + \mu \int |xf'| \, \mathrm{d}x \end{align*}$$ $n$

Mutatis mutandis también puede mostrar que el resultado se mantiene en aritmética de punto fijo.