Gibt es für verrauschte oder fein strukturierte Daten bessere Quadraturen als die Mittelpunktregel?

12

Nur die ersten beiden Abschnitte dieser langen Frage sind wesentlich. Die anderen dienen nur zur Veranschaulichung.

Hintergrund

Fortgeschrittene Quadraturen wie höhergradige zusammengesetzte Newton-Cotes, Gauß-Legendre und Romberg scheinen hauptsächlich für Fälle gedacht zu sein, in denen man die Funktion fein abtasten, aber nicht analytisch integrieren kann. Funktionen mit Strukturen, die kleiner als das Abtastintervall sind (siehe Beispiel in Anhang A) oder Messrauschen können jedoch nicht mit einfachen Ansätzen wie dem Mittelpunkt oder der Trapezregel konkurrieren (siehe Demonstration in Anhang B).

Dies ist etwas intuitiv, da z. B. die zusammengesetzte Simpson-Regel im Wesentlichen ein Viertel der Informationen „verwirft“, indem ein geringeres Gewicht zugewiesen wird. Der einzige Grund, warum solche Quadraturen für ausreichend langweilige Funktionen besser sind, besteht darin, dass der ordnungsgemäße Umgang mit Randeffekten den Effekt verworfener Informationen überwiegt. Aus einer anderen Sicht ist mir intuitiv klar, dass für Funktionen mit feiner Struktur oder Rauschen Proben, die von den Grenzen der Integrationsdomäne entfernt sind, fast gleich weit entfernt sein müssen und fast dasselbe Gewicht haben müssen (für eine hohe Anzahl von Proben) ). Andererseits kann die Quadratur solcher Funktionen von einer besseren Handhabung von Randeffekten profitieren (als bei der Mittelpunktmethode).

Frage

Angenommen, ich möchte verrauschte oder fein strukturierte eindimensionale Daten numerisch integrieren.

Die Anzahl der Abtastpunkte ist fest (aufgrund der hohen Kosten für die Funktionsbewertung), ich kann sie jedoch frei platzieren. Ich (oder die Methode) kann jedoch keine Stichprobenpunkte interaktiv platzieren, dh basierend auf Ergebnissen anderer Stichprobenpunkte. Auch potenzielle Problemregionen kenne ich vorher nicht. So etwas wie Gauß-Legendre (nicht äquidistante Stichprobenpunkte) ist in Ordnung. Adaptive Quadratur ist nicht erforderlich, da interaktiv platzierte Abtastpunkte erforderlich sind.

  • Wurden für einen solchen Fall Methoden vorgeschlagen, die über die Mittelpunktmethode hinausgehen?

  • Oder: Gibt es einen Beweis dafür, dass die Midpoint-Methode unter solchen Bedingungen am besten ist?

  • Allgemeiner: Gibt es Arbeiten zu diesem Problem?

Anhang A: Spezifisches Beispiel einer fein strukturierten Funktion

Ich möchte 1 0 f ( t ) schätzen für: mitund. Eine typische Funktion sieht so aus:01f(t)dtφi[0,2π]logωi[1,1000]

f(t)=ich=1kSünde(ωicht-φich)ωich,
φich[0,2π]Logωich[1,1000]

überlagerte Sinus

Ich habe diese Funktion für die folgenden Eigenschaften gewählt:

  • Es kann für ein Kontrollergebnis analytisch integriert werden.
  • Es hat eine feine Struktur auf einer Ebene, die es unmöglich macht, alles mit der Anzahl der von mir verwendeten Samples zu erfassen ( ).<102
  • Es wird nicht von seiner feinen Struktur dominiert.

Anhang B: Benchmark

Der Vollständigkeit halber ist hier ein Benchmark in 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))

(Ich verwende hier den Median, um den Einfluss von Ausreißern aufgrund von Funktionen mit nur hochfrequenten Inhalten zu verringern. Im Mittel sind die Ergebnisse ähnlich.)

Die Mediane der relativen Integrationsfehler sind:

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

Hinweis: Nach zwei Monaten und einer Prämie ohne Ergebnis habe ich dies auf MathOverflow gepostet .

Wrzlprmft
quelle
Ist dies das Problem, an dem Sie wirklich interessiert sind? In 1D können Sie mit den meisten Methoden wahrscheinlich ziemlich schnell gute Ergebnisse erzielen.
David Ketcheson
"Ich habe eine feste Anzahl von Stichprobenpunkten und kann diese frei platzieren. Ich kann jedoch keine Stichprobenpunkte interaktiv platzieren, dh basierend auf Ergebnissen anderer Stichprobenpunkte." Diese Einschränkung ist mir nicht klar. Darf ich die Knoten dort platzieren, wo ein adaptiver Algorithmus sie platzieren würde, solange ich nur wirklich schlau bin (anstatt den adaptiven Algorithmus tatsächlich zu verwenden)? Wenn ich nicht "wirklich schlau" sein darf, welche Art von Knotenplatzierungen sind dann tatsächlich zulässig?
David Ketcheson
@ DavidKetcheson: Ist dies das Problem, an dem Sie wirklich interessiert sind? - Ja, ich interessiere mich sehr für 1D. - In 1D können Sie mit den meisten Methoden sehr schnell gute Ergebnisse erzielen. - Denken Sie daran, dass die Funktionsbewertung kostspielig sein kann. - Welche Art von Knotenplatzierungen sind dann tatsächlich zulässig? - Ich habe meine Frage bearbeitet, um sie klarer zu machen.
Wrzlprmft
Danke, das hilft. Die Frage erscheint mir immer noch vage. Ich denke, es gibt eine einfache und genauere Frage, die sich besser beantworten lässt. Dazu müssten eine Reihe von Funktionen (die von der zulässigen Anzahl der Quadraturknoten abhängen können) und eine Metrik definiert werden. Dann könnten Sie fragen, ob die Mittelpunktmethode in dieser Metrik für diese Menge von Funktionen optimal ist (wobei vermutlich dieselbe Menge von Knoten für die Quadratur aller Funktionen verwendet werden muss).
David Ketcheson
1
@DavidKetcheson: Hierfür müssten eine Reihe von Funktionen (die von der zulässigen Anzahl der Quadraturknoten abhängen können) und eine Metrik definiert werden. - Da ich zu diesem Thema bisher nichts Nützliches gefunden habe, sehe ich keinen Grund, solche Beschränkungen aufzuerlegen. Vielmehr würde ich mit solchen Einschränkungen das Risiko eingehen, dass ich eine vorhandene Arbeit (oder einen einfachen Beweis) für geringfügig andere Bedingungen oder Annahmen ausschließe. Wenn es Möglichkeiten gibt, das abgebildete Szenario in Definitionen und Ähnlichem zu erfassen, für die es ein Nachschlagewerk oder einen einfachen Beweis gibt, freue ich mich darüber.
Wrzlprmft

Antworten:

1

Zunächst einmal, glaube ich, verstehen Sie das Konzept der adaptiven Quadratur falsch. Adaptive Quadratur impliziert nicht "interaktives Platzieren von Abtastpunkten". Die ganze Idee hinter der adaptiven Quadratur besteht darin, ein Schema zu entwickeln, das eine bestimmte Funktion in einen bestimmten (geschätzten) absoluten oder relativen Fehler mit so wenig Funktionsbewertungen wie möglich integriert.

Eine zweite Bemerkung: Sie schreiben "Die Anzahl der Abtastpunkte ist fest (aufgrund der kostenintensiven Funktionsbewertung), aber ich kann sie frei platzieren". Ich denke, die Idee sollte sein, dass die Anzahl der Abtastpunkte (oder Funktionsbewertungen in Quadratur-Terminologie) so klein wie möglich sein sollte (dh nicht festgelegt).

Welche Idee steckt hinter der adaptiven Quadratur, wie sie beispielsweise in QUADPACK implementiert ist ?

  1. Die Grundzutat ist eine "verschachtelte" Quadraturregel: Dies ist eine Kombination aus zwei Quadraturregeln, bei denen eine höhere Ordnung (oder Genauigkeit) als die andere hat. Warum? Basierend auf dem Unterschied zwischen diesen Regeln kann der Algorithmus den Quadraturfehler schätzen (natürlich verwendet der Algorithmus den genauesten als Referenzergebnis). Beispiele könnten die Trapezregel mit Knoten und Knoten sein. Im Falle von QUADPACK sind die Regeln Gauß-Kronrod-Regeln. Dies sind interpolatorische Quadraturregeln, die eine Gauß-Legendre-Quadraturregel einer bestimmten Ordnung 2 n + 1 N N 2 N - 12n2n+1Nund eine optimale Erweiterung dieser Regel. Dies bedeutet, dass eine höhere Quadraturordnung erhalten werden kann, indem die Gauß-Legendre-Knoten (dh die kostspieligen Funktionsbewertungen) mit unterschiedlichen Gewichten erneut verwendet werden und eine Anzahl zusätzlicher Knoten hinzugefügt wird. Mit anderen Worten, die ursprüngliche Gauß-Legendre-Regel der Ordnung wird alle Polynome des Grades genau integrieren, während die erweiterte Gauß-Kronrod-Regel einige Polynome höherer Ordnung genau integrieren wird. Eine klassische Regel ist die G7K15 (Gauß-Legende 7. Ordnung mit Gauß-Kronrod 15. Ordnung). Die Magie ist, dass die 7 Knoten des Gauss-Legendre eine Teilmenge der 15 Knoten des Gauss-Kronrod sind, also habe ich mit 15 Funktionsbewertungen eine Quadraturbewertung zusammen mit einer Fehlerschätzung!N2N-1

  2. Die nächste Zutat ist eine Strategie zum Teilen und Erobern. Angenommen, Sie lassen diesen G7K15 auf Ihrem Integranden los und stellen einen Quadraturfehler fest, der nach Ihrem Geschmack zu groß ist. QUADPACK unterteilt dann das ursprüngliche Intervall in zwei gleich beabstandete Teilintervalle. Anschließend werden die beiden Teilintegrale mit der Grundregel G7K15 neu ausgewertet. Der Algorithmus hat jetzt eine globale Fehlerschätzung (die hoffentlich niedriger sein sollte als die erste), aber auch zwei lokale Fehlerschätzungen. Es wählt das Intervall mit dem größten Fehler aus und teilt dieses in zwei Teile. Es werden zwei neue Integrale geschätzt und der globale Fehler wird aktualisiert. Und so weiter, bis der globale Fehler unter Ihrem angeforderten Ziel liegt oder die maximale Anzahl von Unterteilungen überschritten wurde.

Daher fordere ich Sie auf, Ihren obigen Code mithilfe der scipy.quadMethode zu aktualisieren . Vielleicht müssen Sie bei einem Integranden mit viel "Feinstruktur" die maximale Anzahl der Unterteilungen erhöhen ( limitOption). Sie können auch mit den epsabsund / oder epsrelParametern spielen.

Wenn Sie jedoch nur experimentelle Daten haben, sehe ich zwei Möglichkeiten.

  1. Wenn Sie die Möglichkeit haben , die Messpunkte, dh Werte auszuwählen , würde ich sie wählen equidistanly und vorzugsweise als eine Potenz von , so dass Sie eine verschachtelte Trapezregel anwenden können (und profitieren Sie von Romberg Extrapolation).2t2
  2. Wenn Sie keine Möglichkeit haben, die Knoten auszuwählen, dh die Messungen zu zufälligen Zeiten erfolgen, ist die beste Option meiner Meinung nach immer noch die Trapezregel.
GertVdE
quelle
Ich denke, Sie verstehen das Konzept der adaptiven Quadratur falsch. - Ihr Beitrag stimmt vollständig mit meinem vorherigen Verständnis der adaptiven Quadratur überein und stimmt eindeutig mit der Definition überein, wie ich die Abtastpunkte interaktiv platziert habe (ob dies nun eine geeignete Phrase ist oder nicht). - Sie schreiben […]. Ich denke, die Idee sollte sein, dass die Anzahl der Stichprobenpunkte […] so gering wie möglich sein sollte (dh nicht festgelegt). - Wenn Sie diesen Luxus haben, sicher, aber experimentelle Einschränkungen können nicht so harmlos sein. Angenommen, Sie müssen etwas gleichzeitig mit einer festen Anzahl teurer Sensoren messen.
Wrzlprmft
Entschuldigen Sie. Ich habe "interaktiv" in Ihrer Frage falsch interpretiert. Nach meinem Verständnis bedeutet "interaktiv", dass der Benutzer eingreift, nicht durch einen Algorithmus. Ich habe einen Absatz in meine Antwort zu experimentellen Daten eingefügt. Ein anderer Ansatz wäre, die Feinstrukturinformationen "herauszufiltern", dh eine Fouriertransformation anzuwenden und Frequenzen höherer Ordnung mit kleinen Amplituden zu entfernen. Wäre das eine Option?
GertVdE
Wenn Sie die Möglichkeit haben, die Messpunkte auszuwählen […] - Äquidistante Punkte sind sowieso das, was ich für Mittelpunkt, einfaches Trapez usw. benötige, also ist dies genau das, was ich in meinem Benchmark gemacht habe. Hier bringt die Romberg-Extrapolation keinen Vorteil.
Wrzlprmft
Ein anderer Ansatz wäre, die Feinstrukturinformationen herauszufiltern […]. Wäre das eine Option? - In meinem Beispiel gehe ich davon aus, dass die Feinstruktur ein Teil dessen ist, was ich messen möchte. Ich habe einfach nicht genügend Samples, um sie vollständig zu erfassen. Was das tatsächliche Rauschen betrifft, gibt es keine technischen Einschränkungen, die mich vom Filtern abhalten. Das Integral über die gesamte Domäne ist jedoch bereits das ultimative Tiefpassfilter. Daher bin ich skeptisch, dass dies verbessert werden kann, ohne dass Rauschen mit bestimmten, harmlosen und bekannten Eigenschaften auftritt.
Wrzlprmft
Ist es wirklich stochastisch? Es muss einige abgeleitete geben, die stochastische Integralnäherungen höherer Ordnung sind.
Chris Rackauckas
0

Ich bin nicht davon überzeugt, dass Ihr Code irgendetwas Grundlegendes über die verschiedenen Quadraturregeln und wie gut sie sich gegen Rauschen und Feinstruktur verhalten, und glaube, dass Sie bei Auswahl verschiedener Feinstrukturen etwas anderes finden würden. Hier ist der Satz:

Keine Quadraturmethode kann einen geringen absoluten oder relativen Fehler gegen eine Funktion mit unbegrenzter Gesamtvariation ergeben. In einem Gleitkommasystem mit Einheitenrundung haben wir die Schätzungμ wobei die Quadratursumme ist, die auf die numerische Implementierung von einwirkt .

|einbfdx-Q.^[f^]||einbfdx-Q.[f]|+μ[4einb|f|dx+einb|xf|dx]
Q.^f^f

Beweis: Die Quadraturknoten seien und die (nicht negativen) Quadraturgewichte seien und bezeichnen ihre Gleitkommanäherungen mit und . Angenommen, erfüllt wobei wobei die Einheitenrundung ist. Dann {xich}ich=0n-1{wich}ich=0n-1w^ichx^ichf^f^(x)=f(x)(1+2δ)|δ|μμ

Q.^[f^]=ich=0n-1w^ichf^(x^ich)=ich=0n-1wich(1+δichw)f(xich+δichxxich)(1+2δichf)(1+δich)ich=0n-1wich[f(xich)+δichxxichf(xich)](1+δichw+2δichf+δich)ich=0n-1wichf(xich)+ich=0n-1δichxwichxichf(xich)+wichf(xich)(δichw+2δichf+δich)
so aus, dass Dies setzt voraus, dass die Summe fehlerfrei berechnet wird. multiplizieren Sie mit , um diese Annahme fallen zu lassen.
|Q.^[f^]-Q.[f]|μich=0n-1wich(|xichf(xich)|+4|f(xich)|)4μ|f|dx+μ|xf|dx
n

Mutatis mutandis können Sie auch zeigen, dass das Ergebnis in Festkomma-Arithmetik gilt.

user14717
quelle
Vielen Dank für Ihre Antwort. Ich habe ein bisschen Probleme, das Szenario zu verstehen, das Sie in Betracht ziehen, und wie es sich auf meine Frage bezieht. Was meinen Sie mit unbegrenzter Gesamtvariation des Gleitkommas? Sofern ich mich nicht sehr irre, sind alle meine Rechenergebnisse (mit Ausnahme des Kontrollfalls mit Romberg und Gauß-Legendre) weit davon entfernt, von Ungenauigkeiten der arithmetischen Implementierung (Gleitkomma oder Festkomma) beeinflusst zu werden. Das Geräusch, das ich betrachte, ist ebenfalls nicht numerischer Natur, sondern experimentell.
Wrzlprmft
@Wrzlprmft: Fließkomma ist das Ergebnis, das ich nachweisen konnte. Ich kann es auch im Festkomma beweisen, was dann anzeigt, dass das Ergebnis für experimentelle Daten gilt. Ich glaube, es gilt für jede Fehlerquelle in den Quadraturknoten. Ich habe zur Verdeutlichung bearbeitet.
user14717
Für experimentelle Daten ist das Ergebnis viel überzeugender, da experimentelle Daten im Allgemeinen nicht differenzierbar sind und daher die Gesamtvariation unendlich ist.
user14717
Es tut mir leid, aber ich folge dir immer noch nicht. Ihr Ergebnis scheint der Fehler zu sein, der bei der numerischen Implementierung der Quadratur gemacht wurde, nicht der Fehler der Quadratur selbst. Das Problem, das ich habe, betrifft das letztere, und insbesondere sehe ich keinen Grund zu der Annahme, dass es sich für nicht manifestieren würde . μ=0
Wrzlprmft
Die Grundidee hierbei ergibt sich aus der Bedingungsnummer der Funktionsbewertung. Ihre Bewertungen sind schlecht konditioniert, da sie laut sind.
user14717