Pour les données bruyantes ou à structure fine, existe-t-il de meilleures quadratures que la règle médiane?

Seules les deux premières sections de cette longue question sont essentielles. Les autres sont juste à titre d'illustration.

Contexte

Des quadratures avancées telles que Newton – Cotes composites de degré supérieur, Gauß – Legendre et Romberg semblent être principalement destinées aux cas où l'on peut finement échantillonner la fonction sans l'intégrer analytiquement. Cependant, pour les fonctions avec des structures plus fines que l'intervalle d'échantillonnage (voir l'exemple A) ou le bruit de mesure, elles ne peuvent pas rivaliser avec des approches simples telles que la règle médiane ou trapézoïdale (voir l'annexe B pour une démonstration).

Ceci est quelque peu intuitif car, par exemple, la règle composite de Simpson «élimine» essentiellement un quart des informations en lui attribuant un poids inférieur. La seule raison pour laquelle ces quadratures sont meilleures pour des fonctions suffisamment ennuyeuses est que la gestion correcte des effets de bordure l'emporte sur l'effet des informations rejetées. D'un autre point de vue, il est intuitivement clair pour moi que pour les fonctions à structure fine ou à bruit, les échantillons éloignés des frontières du domaine d'intégration doivent être quasiment équidistants et avoir presque le même poids (pour un nombre élevé d'échantillons) ). En revanche, la quadrature de ces fonctions peut bénéficier d'une meilleure gestion des effets de frontière (que pour la méthode du point médian).

Question

Supposons que je souhaite intégrer numériquement des données unidimensionnelles bruyantes ou à structure fine.

Le nombre de points d'échantillonnage est fixe (en raison de l'évaluation coûteuse des fonctions), mais je peux les placer librement. Cependant, I (ou la méthode) ne peut pas placer les points d'échantillonnage de manière interactive, c'est-à-dire sur la base des résultats d'autres points d'échantillonnage. Je ne connais pas non plus au préalable les régions à problèmes potentiels. Donc, quelque chose comme Gauß – Legendre (points d'échantillonnage non équidistants) est correct; la quadrature adaptative ne l'est pas car elle nécessite des points d'échantillonnage placés de manière interactive.

• Des méthodes allant au-delà de la méthode du point médian ont-elles été suggérées pour un tel cas?

• Ou: Existe-t-il une preuve que la méthode du point médian est la meilleure dans de telles conditions?

• Plus généralement: existe-t-il des travaux sur ce problème?

Annexe A: Exemple spécifique d'une fonction à structure fine

Je souhaite estimer pour: avec et . Une fonction typique ressemble à ceci:$\int_0^1f(t)\, \mathrm{d}t$φi∈[0,2π]logωi∈[1,1000

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

J'ai choisi cette fonction pour les propriétés suivantes:

• Il peut être intégré analytiquement pour un résultat de contrôle.
• Il a une structure fine à un niveau qui rend impossible de capturer tout cela avec le nombre d'échantillons que j'utilise ( ).$<10^2$
• Il n'est pas dominé par sa structure fine.

Annexe B: Benchmark

Pour être complet, voici une référence 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))

(J'utilise ici la médiane pour réduire l'influence des valeurs aberrantes en raison de fonctions qui n'ont qu'un contenu haute fréquence. Pour la moyenne, les résultats sont similaires.)

Les médianes des erreurs d'intégration relatives sont:

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

Remarque: Après deux mois et une prime sans résultat, j'ai posté ceci sur MathOverflow .

Tout d'abord, je pense que vous comprenez mal le concept de quadrature adaptative. La quadrature adaptative n'implique pas de "placer interactivement des points d'échantillonnage". L'idée derrière la quadrature adaptative est de concevoir un schéma qui intégrera une certaine fonction à une certaine erreur absolue ou relative (estimée) avec aussi peu d'évaluations de fonction que possible.

Une deuxième remarque: vous écrivez "Le nombre de points d'échantillonnage est fixe (car l'évaluation des fonctions est coûteuse), mais je peux les placer librement". Je pense que l'idée devrait être que le nombre de points d'échantillonnage (ou d'évaluations de fonction dans la terminologie en quadrature) soit aussi petit que possible (c'est-à-dire non fixe).

Alors, quelle est l'idée derrière la quadrature adaptative implémentée dans QUADPACK par exemple?

1. L'ingrédient de base est une règle de quadrature «imbriquée»: il s'agit d'une combinaison de deux règles de quadrature où l'une a un ordre (ou une précision) supérieur à l'autre. Pourquoi? En fonction de la différence entre ces règles, l'algorithme peut estimer l'erreur de quadrature (bien sûr, l'algorithme utilisera la plus précise comme résultat de référence). Des exemples pourraient être la règle trapézoïdale avec nœuds et nœuds. Dans le cas de QUADPACK, les règles sont des règles de Gauss-Kronrod. Ce sont des règles de quadrature interpolatoires qui utilisent une règle de quadrature de Gauss-Legendre d'un certain ordre 2 n + 1 N N 2 N - $2^{n}$$2^{n+1}$$N$et une extension optimale de cette règle. Cela signifie qu'un ordre de quadrature supérieur peut être obtenu en réutilisant les nœuds de Gauss-Legendre (c'est-à-dire les évaluations de fonctions coûteuses) avec des poids différents et en ajoutant un certain nombre de nœuds supplémentaires. En d'autres termes, la règle d'origine Gauss-Legendre d'ordre intégrera tous les polynômes de degré exactement tandis que la règle de Gauss-Kronrod étendue intégrera exactement certains polynômes d'ordre supérieur. Une règle classique est le G7K15 (Gauss-Legendre du 7e ordre avec Gauss-Kronrod du 15e ordre). La magie est que les 7 nœuds du Gauss-Legendre sont un sous-ensemble des 15 nœuds du Gauss-Kronrod donc avec 15 évaluations de fonction, j'ai une évaluation en quadrature avec une estimation d'erreur!$N$$2N-1$

2. L'ingrédient suivant est une stratégie de «diviser pour mieux régner». Supposons que vous lâchiez ce G7K15 sur votre intégrande et que vous observiez une erreur de quadrature trop importante selon votre goût. QUADPACK subdivise ensuite l'intervalle d'origine en deux sous-intervalles également espacés. Et puis il réévaluera les deux sous-intégrales en utilisant la règle de base, G7K15. Maintenant, l'algorithme a une estimation d'erreur globale (qui devrait être, espérons-le, inférieure à la première) mais également deux estimations d'erreur locales. Il sélectionne l'intervalle avec l'erreur la plus importante et divise celle-ci en deux. Deux nouvelles intégrales sont estimées et l'erreur globale est mise à jour. Et ainsi de suite jusqu'à ce que l'erreur globale soit inférieure à votre cible demandée ou que le nombre maximum de subdivisions soit dépassé.

Je vous mets donc au défi de mettre à jour votre code ci-dessus en utilisant la scipy.quadméthode. Peut-être que dans le cas d'un integrand avec beaucoup de "structure fine", vous devrez peut-être augmenter le nombre maximum de subdivisions (l' limitoption). Vous pouvez également jouer avec les paramètres epsabset / ou epsrel.

Cependant, si vous ne disposez que de données expérimentales, je vois deux possibilités.

1. Si vous avez la possibilité de sélectionner les points de mesure, c'est-à-dire les valeurs de , je les sélectionnerais équidistamment et de préférence comme une puissance de afin que vous puissiez appliquer une règle trapézoïdale imbriquée (et bénéficier de l'extrapolation de Romberg).$t$$2$
2. Si vous n'avez aucun moyen de choisir les nœuds, c'est-à-dire que les mesures viennent à des moments aléatoires, la meilleure option à mon avis est toujours la règle trapézoïdale.

Je ne suis pas convaincu que votre code démontre quoi que ce soit de fondamental concernant les différentes règles de quadrature et leur efficacité contre le bruit et la structure fine, et je pense qu'il est probable que si vous choisissez différentes structures d'amendes différentes, vous trouverez quelque chose de différent. Voici le théorème:

Aucune méthode en quadrature ne peut donner une faible erreur absolue ou relative par rapport à une fonction à variation totale illimitée. Dans un système à virgule flottante avec arrondi unitaire , nous avons l'estimation$\mu$ où est la somme en quadrature agissant sur l'implémentation numérique 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$

Preuve: que les nœuds de quadrature soient et que les poids de quadrature (non négatifs) soient et désignent leurs approximations à virgule flottante par et . Supposons que satisfait où où est l'arrondi d'unité. alors $\{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*}$$ pour que Cela suppose que la somme est calculée sans erreur; multipliez par pour supprimer cette hypothèse. $$\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, vous pouvez également montrer que le résultat est valable en arithmétique à virgule fixe.