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

12

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í:01f(t)dt

f(t)=i=1ksin(ωitφi)ωi,
φi[0,2π]logωi[1,1000]

senos superpuestos

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 ( ).<102
  • 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 .

Wrzlprmft
fuente
¿Es este el tipo de problema que realmente le interesa? En 1D, probablemente pueda obtener buenos resultados bastante rápido con la mayoría de los métodos.
David Ketcheson el
"Tengo un número fijo de puntos de muestreo y puedo colocarlos libremente. Sin embargo, no puedo colocar puntos de muestreo de manera interactiva, es decir, en función de los resultados de otros puntos de muestreo". Esta restricción no está clara para mí. ¿Se me permite poner los nodos donde los colocaría un algoritmo adaptativo, siempre y cuando sea realmente inteligente (en lugar de usar el algoritmo adaptativo)? Si no se me permite ser "realmente inteligente" al respecto, ¿qué tipo de ubicaciones de nodo se permiten realmente?
David Ketcheson el
@DavidKetcheson: ¿Es este el tipo de problema que realmente le interesa? - Sí, estoy realmente interesado en 1D. - En 1D, probablemente pueda obtener buenos resultados bastante rápido con la mayoría de los métodos. - Recuerde que la evaluación de la función puede ser costosa. - Entonces, ¿qué tipo de ubicaciones de nodos se permiten realmente? - Edité mi pregunta con la esperanza de hacerlo más claro.
Wrzlprmft
Gracias eso ayuda. Para mí, la pregunta todavía parece vaga. Creo que hay una pregunta simple y más precisa que sería más responsable. Requeriría definir un conjunto de funciones (que podrían depender del número permitido de nodos de cuadratura) y una métrica. Luego, podría preguntar si el método del punto medio es óptimo en esa métrica sobre ese conjunto de funciones (donde presumiblemente se debe usar el mismo conjunto de nodos para cuadraturar todas las funciones).
David Ketcheson el
1
@DavidKetcheson: requeriría definir un conjunto de funciones (que podrían depender del número permitido de nodos de cuadratura) y una métrica. - Dado que no he encontrado nada útil sobre este tema hasta el momento, no veo ninguna razón para imponer tales restricciones. Más bien, con tales restricciones, me arriesgaría a excluir algunos trabajos existentes (o pruebas fáciles) por condiciones o suposiciones ligeramente diferentes. Si hay alguna forma de capturar el escenario representado en definiciones y similares para las cuales existe un trabajo de referencia o una prueba fácil, estoy contento con eso.
Wrzlprmft

Respuestas:

1

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 orden2n2n+1Ny 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!N2N1

  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).t2
  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.
GertVdE
fuente
Creo que no entiendes el concepto de cuadratura adaptativa. - Su publicación está completamente de acuerdo con mi comprensión previa de la cuadratura adaptativa y es una coincidencia clara de cómo definí interactivamente colocando los puntos de muestreo (ya sea que sea una frase apropiada o no). - usted escribe [...] Creo que la idea debería ser que el número de puntos de muestreo [...] debería ser lo más pequeño posible (es decir, no fijo). - Si tiene ese lujo, claro, pero las restricciones experimentales pueden no ser tan benignas. Por ejemplo, suponga que tiene que medir algo simultáneamente con un número fijo de sensores caros.
Wrzlprmft
Mis disculpas. Interpreté mal "interactivamente" en su pregunta. En mi opinión, "interactivamente" significa la intervención del usuario, no por un algoritmo. He agregado un párrafo en mi respuesta sobre datos experimentales. Otro enfoque sería "filtrar" la información de estructura fina, es decir, aplicar una transformada de Fourier y eliminar frecuencias de alto orden con pequeñas amplitudes. ¿Sería esa una opción?
GertVdE
Si tiene la oportunidad de seleccionar los puntos de medición […] , de todos modos, los puntos equidistantes son lo que necesito para el punto medio, el trapecio simple, etc., así que esto es exactamente lo que hice en mi punto de referencia. Aquí, la extrapolación de Romberg no ofrece ninguna ventaja.
Wrzlprmft
Otro enfoque sería "filtrar" la información de estructura fina [...] ¿Sería una opción? - En mi ejemplo, supongo que la estructura fina es parte de lo que quiero medir, simplemente no tengo suficientes muestras para capturarla por completo. En cuanto al ruido real, no hay restricciones técnicas que me impidan filtrar. Sin embargo, la integral en todo el dominio ya es el mejor filtro de paso bajo, por lo que soy escéptico de que esto pueda mejorarse sin tener ruido con propiedades específicas, benignas y conocidas.
Wrzlprmft
¿Es realmente estocástico? Debe haber algunos derivados que sean aproximaciones integrales estocásticas de orden superior.
Chris Rackauckas
0

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μ donde es la suma en cuadratura que actúa sobre la implementación numérica de .

|abfdxQ^[f^]||abfdxQ[f]|+μ[4ab|f|dx+ab|xf|dx]
Q^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 {xi}i=0n1{wi}i=0n1w^ix^if^f^(x)=f(x)(1+2δ)|δ|μμ

Q^[f^]=i=0n1w^if^(x^i)=i=0n1wi(1+δiw)f(xi+δixxi)(1+2δif)(1+δi)i=0n1wi[f(xi)+δixxif(xi)](1+δiw+2δif+δi)i=0n1wif(xi)+i=0n1δixwixif(xi)+wif(xi)(δiw+2δif+δi)
para que Esto supone que la suma se calcula sin error; multiplique por para descartar esa suposición.
|Q^[f^]Q[f]|μi=0n1wi(|xif(xi)|+4|f(xi)|)4μ|f|dx+μ|xf|dx
n

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

usuario14717
fuente
Gracias por tu respuesta. Tengo algunos problemas para comprender el escenario que está considerando y cómo se relaciona con mi pregunta. ¿Qué quiere decir con variación total ilimitada en coma flotante? A menos que esté muy equivocado, todos mis resultados computacionales (excepto el caso de control con Romberg y Gauß – Legendre) están lejos de estar influenciados por imprecisiones de la implementación aritmética (punto flotante o punto fijo). El ruido que estoy considerando tampoco es de naturaleza numérica, sino experimental.
Wrzlprmft
@Wrzlprmft: El punto flotante es el resultado que pude probar. También puedo probarlo en un punto fijo, que luego indica que el resultado se cumple para datos experimentales. Creo que es cierto para cualquier fuente de error en los nodos de cuadratura. He editado para aclarar.
user14717
Para los datos experimentales, el resultado es mucho más convincente porque, en general, los datos experimentales no son diferenciables y, por lo tanto, la variación total es infinita.
user14717
Lo siento, pero todavía no puedo seguirte. Su resultado parece ser sobre el error cometido al implementar numéricamente la cuadratura, no sobre el error de la propia cuadratura. El problema que tengo es sobre esto último y, en particular, no veo ninguna razón para creer que no se manifestaría para . μ=0
Wrzlprmft
La idea principal aquí proviene del número de condición de la evaluación de la función. Sus evaluaciones están mal acondicionadas ya que son ruidosas.
user14717