Czy w przypadku danych zaszumionych lub o drobnej strukturze są lepsze kwadratury niż reguła punktu środkowego?

12

Niezbędne są tylko dwie pierwsze części tego długiego pytania. Pozostałe służą wyłącznie ilustracji.

tło

Zaawansowane kwadratury, takie jak kompozyt wyższego stopnia Newtona-Cotesa, Gaußa-Legendre'a i Romberga, wydają się być głównie przeznaczone do przypadków, w których można dokładnie próbkować funkcję, ale nie integrować analitycznie. Jednak w przypadku funkcji o strukturach większych niż interwał próbkowania (patrz przykład A dodatek) lub szumu pomiarowego nie mogą one konkurować z prostymi podejściami, takimi jak punkt środkowy lub reguła trapezowa (demonstracja znajduje się w dodatku B).

Jest to nieco intuicyjne, ponieważ np. Złożona reguła Simpsona zasadniczo „odrzuca” jedną czwartą informacji, przypisując jej niższą wagę. Jedynym powodem, dla którego takie kwadratury są lepsze dla wystarczająco nudnych funkcji, jest to, że prawidłowe zarządzanie efektami granicznymi przeważa nad efektem odrzucenia informacji. Z innego punktu widzenia intuicyjnie jest dla mnie jasne, że w przypadku funkcji o drobnej strukturze lub szumie próbki, które są oddalone od granic domeny integracji, muszą być prawie w równej odległości i mieć prawie taką samą wagę (w przypadku dużej liczby próbek ). Z drugiej strony kwadratura takich funkcji może skorzystać na lepszej obsłudze efektów granicznych (niż w metodzie punktu środkowego).

Pytanie

Załóżmy, że chcę zintegrować numerycznie jednowymiarowe dane z zaszumieniem lub drobną strukturą.

Liczba punktów próbkowania jest stała (ze względu na kosztowną ocenę funkcji), ale mogę je dowolnie umieszczać. Jednak ja (lub metoda) nie mogę umieszczać punktów próbkowania interaktywnie, tj. W oparciu o wyniki z innych punktów próbkowania. Nie znam też wcześniej potencjalnych regionów problemowych. Tak więc coś w stylu Gauß – Legendre (nierównoodległe punkty próbkowania) jest w porządku; kwadratura adaptacyjna nie jest, ponieważ wymaga interaktywnie rozmieszczonych punktów próbkowania.

  • Czy w takim przypadku zaproponowano metody wykraczające poza metodę punktu środkowego?

  • Lub: Czy istnieje dowód, że metoda punktu środkowego jest najlepsza w takich warunkach?

  • Mówiąc bardziej ogólnie: czy są jakieś prace nad tym problemem?

Dodatek A: Konkretny przykład funkcji o drobnej strukturze

Chcę oszacować dla: z i . Typowa funkcja wygląda następująco:01f(t)dtφi[0,2π]logωi[1,1000]

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

nałożone sinus

Wybrałem tę funkcję dla następujących właściwości:

  • Można go zintegrować analitycznie w celu uzyskania wyniku kontroli.
  • Ma drobną strukturę na poziomie, który uniemożliwia uchwycenie wszystkiego przy użyciu liczby próbek, których używam ( ).<102
  • Nie jest zdominowany przez swoją delikatną strukturę.

Załącznik B: Benchmark

Dla kompletności, oto test porównawczy w Pythonie:

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))

(Używam tutaj mediany, aby zmniejszyć wpływ wartości odstających z powodu funkcji, które mają tylko zawartość wysokich częstotliwości. Dla średniej wyniki są podobne.)

Mediany względnych błędów integracji to:

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

Uwaga: po dwóch miesiącach i jednej nagrodie bez rezultatu opublikowałem to na MathOverflow .

Wrzlprmft
źródło
Czy jest to problem, który naprawdę Cię interesuje? W 1D prawdopodobnie możesz uzyskać dobre wyniki dość szybko za pomocą większości metod.
David Ketcheson
„Mam określoną liczbę punktów próbkowania i mogę je dowolnie umieszczać. Nie mogę jednak umieszczać punktów próbkowania interaktywnie, tj. W oparciu o wyniki z innych punktów próbkowania”. To ograniczenie nie jest dla mnie jasne. Czy mogę umieszczać węzły tam, gdzie umieściłby je algorytm adaptacyjny, o ile jestem po prostu naprawdę inteligentny (zamiast faktycznie używać algorytmu adaptacyjnego)? Jeśli nie wolno mi być „naprawdę mądrym” w tym zakresie, to jakie umiejscowienie węzłów jest w rzeczywistości dozwolone?
David Ketcheson
@DavidKetcheson: Czy jest to problem, który naprawdę Cię interesuje? - Tak, naprawdę interesuje mnie 1D. - W 1D prawdopodobnie możesz uzyskać dobre wyniki dość szybko za pomocą większości metod. - Pamiętaj, że ocena funkcji może być kosztowna. - w takim razie jakiego rodzaju umiejscowienie węzłów jest faktycznie dozwolone? - Zredagowałem moje pytanie, mając nadzieję, że wyjaśnię je bardziej.
Wrzlprmft
Dzięki, że pomaga. Dla mnie pytanie wciąż wydaje się niejasne. Myślę, że istnieje proste i bardziej precyzyjne pytanie, na które można by bardziej odpowiedzieć. Wymagałoby to zdefiniowania zestawu funkcji (które mogą zależeć od dozwolonej liczby węzłów kwadraturowych) i metryki. Następnie możesz zapytać, czy metoda punktu środkowego jest optymalna w tej metodzie w porównaniu z tym zestawem funkcji (gdzie prawdopodobnie ten sam zestaw węzłów musi być użyty do kwadratury wszystkich funkcji).
David Ketcheson
1
@DavidKetcheson: Wymagałoby to zdefiniowania zestawu funkcji (które mogą zależeć od dozwolonej liczby węzłów kwadraturowych) i metryki. - Biorąc pod uwagę, że jak dotąd nie znalazłem nic przydatnego w tym temacie, nie widzę powodu, by nakładać takie ograniczenia. Raczej przy takich ograniczeniach zaryzykowałbym wykluczenie niektórych istniejących prac (lub łatwych dowodów) dla nieco innych warunków lub założeń. Jeśli istnieją sposoby uchwycenia przedstawionego scenariusza w definicjach i podobnych, dla których istnieje praca referencyjna lub łatwy dowód, cieszę się z tego.
Wrzlprmft

Odpowiedzi:

1

Po pierwsze, myślę, że źle rozumiesz pojęcie kwadratury adaptacyjnej. Kwadratura adaptacyjna nie oznacza „interaktywnego umieszczania punktów próbki”. Cała idea kwadratury adaptacyjnej polega na opracowaniu schematu, który zintegruje określoną funkcję z określonym (szacowanym) błędem bezwzględnym lub względnym przy możliwie najmniejszej ocenie funkcji.

Druga uwaga: piszesz „Liczba punktów próbkowania jest stała (ze względu na kosztowną ocenę funkcji), ale mogę je dowolnie umieszczać”. Myślę, że powinna istnieć idea, aby liczba punktów próbkowania (lub ocen funkcji w terminologii kwadraturowej) była jak najmniejsza (tj. Nie ustalona).

Jaki jest więc pomysł na kwadraturę adaptacyjną zaimplementowaną na przykład w QUADPACK ?

  1. Podstawowym składnikiem jest „zagnieżdżona” zasada kwadraturowa: jest to kombinacja dwóch reguł kwadraturowych, w których jedna ma wyższą kolejność (lub dokładność) jak druga. Dlaczego? Na podstawie różnicy między tymi regułami algorytm może oszacować błąd kwadraturowy (oczywiście algorytm użyje najbardziej dokładnego jako wyniku odniesienia). Przykładami mogą być reguła trapezowa z węzłami i węzłami. W przypadku QUADPACK są to reguły Gaussa-Kronroda. Są to interpolacyjne reguły kwadraturowe, które wykorzystują regułę kwadraturową Gaussa-Legendre'a określonego rzędu 2 n + 1 N N 2 N - 12n2n+1Ni optymalne rozszerzenie tej zasady. Oznacza to, że można uzyskać wyższą kolejność kwadraturową, ponownie wykorzystując węzły Gaussa-Legendre'a (tj. Kosztowne oceny funkcji) o różnych wagach i dodając szereg dodatkowych węzłów. Innymi słowy, oryginalna reguła Gaussa-Legendre'a rzędu zintegruje dokładnie wszystkie wielomiany stopnia dokładnie, podczas gdy rozszerzona reguła Gaussa-Kronroda dokładnie zintegruje jakiś wielomian wyższego rzędu. Klasyczną zasadą jest G7K15 (Gauss-Legendre 7. rzędu z Gauss-Kronrod 15. rzędu). Magia polega na tym, że 7 węzłów Gaussa-Legendre'a jest podzbiorem 15 węzłów Gaussa-Kronroda, więc z 15 ocenami funkcji mam ocenę kwadratury wraz z oszacowaniem błędu!N2N1

  2. Kolejnym składnikiem jest strategia „dziel i rządź”. Załóżmy, że puściłeś G7K15 na swoim integrandzie i zaobserwowałeś błąd kwadratury, który według twojego gustu jest zbyt duży. QUADPACK podzieli pierwotny interwał na dwie równomiernie rozmieszczone podinterwały. A następnie dokona ponownej oceny dwóch subintegrałów przy użyciu podstawowej zasady, G7K15. Teraz algorytm ma globalne oszacowanie błędów (które powinno być, być może, niższe niż pierwsze), ale także dwa lokalne oszacowania błędów. Wybiera przedział z największym błędem i dzieli go na dwa. Oszacowano dwie nowe całki i zaktualizowano błąd globalny. I tak dalej, dopóki błąd globalny nie spadnie poniżej żądanego celu lub nie zostanie przekroczona maksymalna liczba podziałów.

Wzywam więc do aktualizacji powyższego kodu za pomocą tej scipy.quadmetody. Być może w przypadku integrandu z dużą liczbą „drobnych struktur” może być konieczne zwiększenie maksymalnej liczby podziałów ( limitopcja). Możesz także grać z parametrami epsabsi / lub epsrel.

Jeśli jednak masz tylko dane eksperymentalne, widzę dwie możliwości.

  1. Jeśli masz możliwość wybrania punktów pomiarowych, tj. Wartości , wybrałbym je jednakowo i najlepiej jako potęgę , abyś mógł zastosować zagnieżdżoną regułę trapezoidalną (i czerpać zyski z ekstrapolacji Romberga).2t2
  2. Jeśli nie masz możliwości wyboru węzłów, tzn. Pomiary są wykonywane losowo, najlepszą moim zdaniem nadal jest reguła trapezowa.
GertVdE
źródło
Myślę, że źle rozumiesz pojęcie kwadratury adaptacyjnej. - Twój post całkowicie zgadza się z moim wcześniejszym rozumieniem kwadratury adaptacyjnej i jest to wyraźne dopasowanie do tego, jak zdefiniowałem interaktywnie umieszczanie punktów próbkowania (czy to jest odpowiednie zdanie, czy nie). - piszesz […]. Myślę, że należy polegać na tym, aby liczba punktów próbkowania […] była jak najmniejsza (tj. Nie ustalona). - Oczywiście, jeśli masz ten luksus, ale eksperymentalne ograniczenia mogą nie być tak łagodne. Załóżmy na przykład, że musisz mierzyć coś jednocześnie za pomocą stałej liczby drogich czujników.
Wrzlprmft
Przepraszam. Źle zinterpretowałem „interaktywnie” w twoim pytaniu. W moim rozumieniu „interaktywnie” oznacza interwencję użytkownika, a nie algorytm. W odpowiedzi dodałem akapit dotyczący danych eksperymentalnych. Innym podejściem byłoby „odfiltrowanie” drobnej struktury informacji, tj. Zastosowanie transformacji Fouriera i usunięcie wysokich częstotliwości z małymi amplitudami. Czy to byłaby opcja?
GertVdE
Jeśli masz możliwość wyboru punktów pomiarowych […] - W każdym razie potrzebuję punktów równych dla punktu środkowego, zwykłego trapezu itp., Więc dokładnie to zrobiłem w swoim teście. Ekstrapolacja Romberga nie przynosi tutaj żadnych korzyści.
Wrzlprmft,
Innym podejściem byłoby „odfiltrowanie” drobnej struktury informacji […] Czy byłaby to opcja? - W moim przykładzie zakładam, że drobna struktura jest częścią tego, co chcę zmierzyć, po prostu nie mam wystarczającej liczby próbek, aby ją całkowicie uchwycić. Jeśli chodzi o rzeczywisty hałas, nie ma technicznych ograniczeń, które powstrzymywałyby mnie przed filtrowaniem. Jednak całka w całej domenie jest już ostatecznym filtrem dolnoprzepustowym, więc jestem sceptyczny, że można to poprawić bez szumu o określonych, łagodnych i znanych właściwościach.
Wrzlprmft,
Czy to naprawdę stochastyczne? Muszą być pewne pochodne, które są przybliżeniami całki stochastycznej wyższego rzędu.
Chris Rackauckas,
0

Nie jestem przekonany, że twój kod pokazuje coś fundamentalnego w różnych regułach kwadratury i jak dobrze radzą sobie z hałasem i drobną strukturą, i wierzę, że jeśli wybierzesz różne struktury drobnych kar, znajdziesz coś innego. Oto twierdzenie:

Żadna metoda kwadraturowa nie może dać niskiego błędu bezwzględnego lub względnego względem funkcji z nieograniczoną zmiennością całkowitą. W systemie zmiennoprzecinkowym z zaokrągleniem jednostek mamy oszacowanieμ gdzie jest sumą kwadraturową działającą na implementacji numerycznej z .

|abfdxQ^[f^]||abfdxQ[f]|+μ[4ab|f|dx+ab|xf|dx]
PQ^f ff^f

Dowód: niech węzły kwadraturowe będą mieć a (nieujemne) wagi kwadraturowe będą wynosić i określają przybliżone liczby zmiennoprzecinkowe przez i . Załóżmy, że spełnia gdzie gdzie jest jednostkowym zaokrągleniem. Następnie {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)
, aby Zakłada się, że suma jest obliczana bez błędów; pomnóż przez aby odrzucić to założenie.
|Q^[f^]Q[f]|μi=0n1wi(|xif(xi)|+4|f(xi)|)4μ|f|dx+μ|xf|dx
n

Mutatis mutandis można również pokazać, że wynik jest zachowany w arytmetyki punktu stałego.

użytkownik14717
źródło
Dziękuję za odpowiedź. Mam trochę problemów ze zrozumieniem rozważanego scenariusza i jego związku z moim pytaniem. Co rozumiesz przez nieograniczoną całkowitą zmienność zmiennoprzecinkową? O ile się nie mylę, wszystkie moje wyniki obliczeniowe (z wyjątkiem przypadku kontrolnego z Rombergiem i Gauß-Legendre) są dalekie od wpływu niedokładności implementacji arytmetycznej (zmiennoprzecinkowej lub stałej). Hałas, który rozważam, nie ma również charakteru liczbowego, ale eksperymentalny.
Wrzlprmft
@Wrzlprmft: Zmienny punkt to wynik, który udało mi się udowodnić. Mogę również udowodnić to w ustalonym punkcie, co następnie wskazuje, że wynik dotyczy danych eksperymentalnych. Uważam, że jest to prawdą w przypadku każdego źródła błędu w węzłach kwadraturowych. Zredagowałem, aby wyjaśnić.
user14717
W przypadku danych eksperymentalnych wynik jest znacznie bardziej przekonujący, ponieważ ogólnie dane eksperymentalne nie są zróżnicowane, a zatem całkowita zmienność jest nieskończona.
user14717
Przykro mi, ale nadal nie mogę Cię śledzić. Twój wynik wydaje się dotyczyć błędu popełnionego podczas numerycznej implementacji kwadratury, a nie błędu samego kwadratury. Problem, który mam, dotyczy tego ostatniego, a w szczególności nie widzę powodu, by sądzić, że nie przejawiałoby się to dla . μ=0
Wrzlprmft
Główna idea tutaj pochodzi z warunku oceny funkcji. Twoje oceny są źle uwarunkowane, ponieważ są głośne.
user14717