Zagadnienia

15. Metody quasi-Monte Carlo

15.1. Co to są metody quasi-Monte Carlo?

Metody Monte Carlo potrafią pokonać przekleństwo wymiaru, mają jednak również swoje wady. Najważniejsze z nich to:

  • (i) niepewność wyniku (który ma charakter probablistyczny),

  • (ii) stosunkowo wolna zbieżność (praktycznie N-1/2),

  • (iii) konieczność stosowania (niekiedy skomplikowanych) generatorów liczb losowych.

Można powiedzieć, że skoro używamy z powodzeniem generatorów liczb pseudo-losowych, to implementacje Monte Carlo są w istocie deterministyczne. Dlatego, wybierając punkty deterministycznie, ale tak, aby w jakiś sposób ,,udawały” i ,,uśredniały” wybór losowy, powinniśmy dostać metodę deterministyczną o zbieżności co najmniej N-1/2. A przy okazji usunęlibyśmy dwie z wymienionych wad Monte Carlo.

Takie intuicyjne myślenie wydaje się nie mieć racjonalnych podstaw, bo przecież metody deterministyczne podlegają przekleństwu wymiaru. Pamiętajmy jednak, że twierdzenie 13.4 o przekleństwie zachodzi w klasie Fr⁢D funkcji różniczkowalnych r razy po każdej zmiennej. Zwróciliśmy na to uwagę już na początku rozdziału 14. Dlatego zasadne jest poszukiwanie pozytywnych rozwiązań dla funkcji o innej regularności.

Quasi-Monte Carlo jest deterministycznym odpowiednikiem klasycznej metody Monte Carlo dla aproksymacji całek na kostkach,

Id⁢f⁢:=⁢∫Df⁢x→⁢⁢d⁢x→,⁢D=0,1d.
Definicja 15.1

Metoda quasi-Monte Carlo polega na przybliżeniu Id⁢f średnią arytmetyczną,

QMCd,N⁢f⁢=⁢QMCd,N⁢t→1,…,t→N⁢:=⁢1n⋅∑j=1Nf⁢t→j,

gdzie t→1,t→2,…,t→N są pewnymi szczególnymi punktami w D wybranymi w sposób deterministyczny.

Przez długi czas od swojego powstania uważano, że metody quasi-Monte Carlo są efektywne jedynie dla całek o niskich wymiarach. Dopiero pod koniec lat 90-tych ubiegłego wieku zauważono, że dają istotnie lepsze wyniki niż Monte Carlo w obliczaniu wartości niektórych instrumentów finansowych. Obecnie quasi-Monte Carlo jest powszchnie uznaną i bardzo popularną metodą, której znaczenie trudno przecenić mimo, że dotychczas nie udało się znaleźć pełnego teoretycznego wytłumaczenia jej efektywności.

15.2. Dyskrepancja

Rozważania na temat quasi-Monte Carlo zaczniemy od zdefiniowania pojęcia dyskrepancji, które odgrywa fundamentalną rolę w analizie efektywności tych metod.

Dla x→=x1,…,xdT∈D niech

0→,x→⁢:=×0,x1×⋯×0,xd

oznacza d wymiarowy prostokąt w D, ,,zakotwiczony” w zerze. Zakładamy przy tym, że 0→,0→=∅.

Definicja 15.2

Dyskrepancją (z gwiazdką) danego zbioru N punktów t→j∈0,1d, 1≤j≤N, nazywamy wielkość

DISCd*⁢t→1,…,t→N⁢=⁢supx→∈D⁡⁢⁢DISCd⁢x→;t→1,…,t→N⁢,

gdzie

DISCd⁢x→;t→1,…,t→N⁢=⁢x1⁢…⁢xd-#⁢j⁢:⁢t→j∈0→,x→N,

a #⁢A oznacza liczbę elementów zbioru A.

Ponieważ x1⁢⋯⁢xd jest d-wymiarową objętością prostokąta 0→,x→, dyskrepancja pokazuje jak dobrze danych N punktów aproksymuje objętości tych prostokątów. Równoważnie, oznaczając przez 10→,x→ funkcję charakterystyczną prostokąta 0→,x→ mamy

∫D10→,x→⁢t⁢⁢d⁢t=x1⁢⋯⁢xd⁢⁢oraz⁢⁢1N⁢∑j=1N10→,x→⁢t→j=#⁢j⁢:⁢t→j∈0→,x→N.

Stąd dyskrepancję można również interpretować jako maksymalny błąd aproksymacji funkcji charakterystycznych prostokątów przez odpowiedni algorytm quasi-Monte Carlo.

Pojęcie dyskrepancji zilustrujemy najpierw na przykładzie jednowymiarowym d=1. Nietrudno pokazać, że dla dowolnych tj

DISC1*⁢t1,…,tN⁢≥⁢12⁢N.(15.1)

Rzeczywiście, DISC1⁢x;t1,…,tN jako funkcja x ma na każdym z przedziałów 0,t1, tj-1,tj, 2≤j≤N, pochodną równą 1, oraz przyjmuje zero dla x=0. Zatem, jeśli dyskrepancja byłaby mniejsza od 1/2⁢N to mielibyśmy

t1<12⁢N,⁢tj-tj-1<1N,⁢2≤j≤N,

a stąd

tN⁢=⁢t1+∑j=2Ntj-tj-1⁢<⁢12⁢N+N-1N⁢= 1-12⁢N.

Otrzymujemy sprzeczność, bo dla tN<x<1-1/2⁢N

DISC1⁢x;t1,…,tN⁢<⁢1-12⁢N-1⁢=-12⁢N.

Z dowodu wynika, że równość w (15.1) zachodzi jedynie dla równomiernego rozmieszczenia punktów,

tj*⁢=⁢j-1/2N,⁢1≤j≤N.

W tym przypadku algorytm QMC1,N redukuje się do zasady punktu środkowego.

Z punktu widzenia praktycznych obliczeń dobrze byłoby mieć ciąg nieskończony

t1,t2,…,tn,…⁢⊂0,1

i konstruować kolejne aproksymacje używając N początkowych wyrazów tego ciągu. Ciekawe, że wtedy zbieżność N-1 nie może być zachowana. Dokładniej, można pokazać istnienie c>0 takiego, że dla każdego ciągu nierówność

DISC1*⁢t1,…,tN⁢≥⁢c⁢⁢ln⁡NN

zachodzi dla nieskończenie wielu N.

Analiza dyskrepancji w wymiarach d≥2 just dużo bardziej skomplikowana. Na razie ograniczymy się do zauważenia, że siatka równomierna jest fatalnym wyborem. Dokładniej, niech

t→j1,…,jd⁢=⁢tj1*,…,tjd*T,⁢1≤ji≤n,1≤i≤d,

będzie równomiernym rozkładem N=nd punktów w D. Rozpatrzenie prostokąta

0,12⁢n×0,1×⋯×0,1︸d-1

wystarcza, aby się przekonać, że wtedy dyskrepancja wynosi co najmniej 12⁢n=12⁢N-1/d.

15.3. Błąd quasi-Monte Carlo

15.3.1. Formuła Zaremby

Jeśli f:0,1→R jest funkcją różniczkowalną to dla każdego x mamy

f⁢x⁢=⁢f⁢1-∫x1f′⁢t⁢⁢d⁢t⁢=⁢f⁢1-∫011x,1⁢t⁢⁢f′⁢t⁢⁢d⁢t.(15.2)

Wzór ten uogólnimy na przypadek funkcji wielu zmiennych w następujący sposób.

Najpierw wprowadzimy kilka niezbędnych oznaczeń. Dla podzbioru indeksów U,

∅⁢≠⁢U⁢⊆⁢1,2,…,d,

niech DU=0,1U, gdzie U jest liczbą elementów w U. Niech dalej x→U∈DU będzie wektorem powstałym z wektora x→=x1,x2,…,xdT∈D poprzez usunięcie współrzędnych xj z j∉U, a x→U;1∈D wektorem, którego j-ta współrzędna wynosi xj dla j∈U oraz wynosi 1 dla j∉U. Na przykład, jeśli d=5 i U=1,4 to dla x→=x1,x2,x3,x4,x5T mamy x→U=x1,x4T i x→U;1=x1,1,1,x4,1T.

W końcu, niech

fU′⁢=⁢∂U⁡f∏j∈U∂⁡uj

będzie skrótowym zapisem odpowiedniej pochodnej mieszanej funkcji f.

Lemat 15.1

Jeśli funkcja f:D→R ma ciągłe pochodne cząstkowe mieszane fU′ dla wszystkich ∅≠U⊆1,2,…,d to

f⁢x→⁢=⁢f⁢1→⁢+⁢∑U-1U⁢∫DU1x→U,1→U⁢z→U⁢⁢fU′⁢z→U;1⁢⁢d⁢z→U,⁢x→∈D(15.3)

(1→=1,1,…,1∈Rd).

Dowód

Dowód przeprowadzimy przez indukcję względem d. Dla d=1 równość (15.3) jest równoważna (15.2). Niech więc d≥2. Stosując założenie indukcyjne do f z ustaloną ostatnią współrzędną xd mamy

f⁢x→⁢=⁢f⁢x→d;1⁢+⁢∑V-1V⁢∫DV1x→V,1→V⁢z→V⁢⁢fV′⁢z→V,xd;1⁢⁢d⁢z→V,

gdzie sumowanie jest po wszystkich ∅≠V⊆1,2,…,d-1. Stosując dalej wzór (15.2) ze względu na xd odpowiednio do f⁢x→d;1 i fV′⁢z→V,xd;1 dostajemy

f⁢x→=f⁢1→⁢-⁢∫Dd1x→d,1→d⁢z→d⁢⁢fd′⁢z→d;1⁢⁢d⁢z→d
+-1V⁢∫DV1x→V,1→V⁢z→V⁢⁢fV′⁢z→V;1⁢⁢d⁢z→V
+∑V∪d-1V∪d⁢∫DV∪d1x→V∪d,1→V∪d⁢z→V∪d⁢⁢fV∪d′⁢z→V∪d;1⁢⁢d⁢z→V∪d.

Teraz wystarczy porównać otrzymaną formułę do prawej strony (15.3). Drugi składnik odpowiada U=d, trzeci składnik podzbiorom U≠∅ do których nie należy d, a czwarty podzbiorom U takim, że d∈U i U≥2.

∎

Zauważmy, że rozwijając f⁢x→ i f⁢t→j zgodnie ze wzorem (15.3) mamy

∫Df⁢x→⁢⁢d⁢x→⁢=⁢f⁢1→+∑U-1U⁢∫DU∫D1x→U,1→U⁢z→U⁢⁢d⁢x→⁢⁢fU′⁢z→U;1⁢⁢d⁢z→U,
1N⁢⁢∑j=1Nf⁢t→j⁢=⁢f⁢1→+∑U-1U⁢∫DU1N⁢∑j=1N1t→jU,1→U⁢z→U⁢⁢fU′⁢z→U;1⁢⁢d⁢z→U.

Ponieważ wartość funkcji charakterystycznej odcinka a→,b→ w punkcie c→ jest równa wartości funkcji charakterystycznej odcinka 0→,c→ w punkcie a→ to

∫D1x→U,1→U⁢z→U⁢⁢d⁢x→⁢=⁢∫DU10→U,z→U⁢x→U⁢⁢d⁢x→U⁢=⁢∏j∈Uzj,
1N⁢⁢∑j=1N1t→jU,1→U⁢z→U⁢=⁢1N⁢⁢#⁢j⁢:⁢t→jU∈0→U,z→U⁢=⁢1N⁢⁢#⁢j⁢:⁢t→j∈0→,z→U;1.

Stąd otrzymujemy następującą formułę Zaremby na błąd quasi-Monte Carlo:

Id⁢f-QMCd,N⁢f⁢=⁢∑U-1U⁢∫DUDISCd⁢z→U;1;t→1,…,t→N⋅fU′⁢z→U;1⁢⁢d⁢z→U.(15.4)

15.3.2. Nierówność Koksmy-Hlawki

Oznaczmy przez Vd⁢D klasę funkcji f:D→R, których pochodne mieszane fU′ istnieją i są ciągłe dla wszystkich ∅≠U⊆1,…,d.

Definicja 15.3

Wahaniem (w sensie Hardy-Krause) funkcji f∈Vd⁢D nazywamy wielkość

Vd⁢f⁢:=⁢∑U∫DUfU′⁢z→U;1⁢⁢d⁢z→U.

Zauważmy, że dla d=1,

V1⁢f⁢=⁢∫01f′⁢x⁢⁢d⁢x⁢=⁢sup⁡⁢∑j=1kf⁢zj-f⁢zj-1⁢:⁢k≥1,0=z0<z1<…<zk=1⁢

jest wahaniem funkcji w klasycznym sensie. Następujące bardzo ważne oszacowanie błędu metody quasi-Monte Carlo nosi nazwę nierówności Koksmy-Hlawki.

Twierdzenie 15.1

Jeśli f∈Vd⁢D to błąd metody quasi-Monte Carlo

QMCd,N⁢f=1N⁢∑j=1Nf⁢t→j

dla aproksymacji całki Id=∫Df⁢x→⁢⁢d⁢x→ szacuje się przez

Id⁢f-QMCd,N⁢f⁢≤⁢DISCd*⁢t→1,…,t→N⋅Vd⁢f.
Dowód

Nierówność wynika natychmiast z formuły Zaremby (15.4), bowiem

Id⁢f-QMCd,N⁢f≤∑U∫DUDISCd⁢z→U;1;t→1,…,t→N⋅fU′⁢z→U;1⁢⁢d⁢z→U
≤DISCd*⁢t→1,…,t→N⋅∑U∫DUfU′⁢z→U;1⁢⁢d⁢z→U.
∎

W nierówności Koksmy-Hlawki czynnik błędu Vd⁢f zależny jedynie od funkcji jest oddzielony od czynnika błędu DISCd*⁢t→1,…,t→N zależnego od wyboru punktów. O ile nie mamy wpływu na wahanie funkcji, możemy starać się wybrać punkty t→j tak, aby zminimalizować ich dyskrepancję. Zasadnicze pytanie brzmi: jak mała może być dyskrepacja? W szczególności, czy można wybrać nieskończony ciąg punktów tak, że oparte na nim algorytmy quasi-Monte Carlo pokonują przekleństwo wymiaru w klasie Vd⁢D?

Na miejscu jest teraz uwaga, że dzięki obecności pochodnych mieszanych rozumowanie analogiczne do tego z dowodu twierdzenia 13.4 prowadzi w klasie funkcji f∈Vd⁢D z V⁢f≤1 jedynie do oszacowania z dołu błędu przez c⁢N-1 (a nie c⁢N-r/d).

Okazuje się, że pełne odpowiedzi na zadane pytania nie są znane. Najlepsze ciągi nieskończone t→1⁢,t→2⁢,…,t→n⁢,… spełniają nierówność

DISCd*⁢t→1⁢,…,t→N⁢⁢≤⁢Cd⁢⁢lnd⁡NN,⁢N=1,2,3,…,

gdzie Cd>0 nie zależy od N.

Wydaje się więc, że w klasie Vd⁢D metody quasi-Monte Carlo dają istotnie lepsze wyniki od Monte Carlo, bo błąd nie tylko jest deterministyczny, ale też dla f∈Vd⁢D zbiega do zera dużo szybciej, tzn. błąd jest rzędu N-1+ϵ dla dowolnego ϵ>0 w przypadku quasi-Monte Carlo, oraz N-1/2 w przypadku Monte Carlo. Nie do końca jest to prawdą. Zauważmy bowiem, że w praktycznych obliczeniach wnioski asymptotyczne nie mają zastosowania, gdy wymiar d jest bardzo duży. Wtedy czynnik Cd⁢lnd⁡N może mieć istotne znaczenie i wręcz sprawiać, że nierówność Koksmy-Hlawki staje się bezużyteczna. Dodatkowo, dobre oszacowanie Cd jest wyjątkowe trudne.

A jednak metody quasi-Monte Carlo są z powodzeniem stosowane w praktyce obliczeniowej nawet dla dużych wymiarów d. Istnieje wiele hipotez tworzonych w celu wyjaśnienia tego pozornego paradoksu. Jedna z najbardziej popularnych i już dość dobrze uzasadnionych teoretycznie mówi, że w praktyce mamy co prawda do czynienia z funkcjami bardzo wielu zmiennych, ale istotnych jest jedynie kilka zmiennych albo grupy kilku zmiennych. Matematycznie oznacza to, że w odpowiadających tym założeniom przestrzeniach zachodzą dużo mocniejsze odpowiedniki nierówności Koksmy-Hlawki.

15.4. Ciągi o niskiej dyskrepancji

Definicja 15.4

Ciąg nieskończony t→1,t→2,t→3,… nazywamy ciągem o niskiej dyskrepancji jeśli istnieje Cd>0 taka, że dla wszystkich N

DISCd*⁢t→1,…,t→N⁢≤⁢Cd⁢⁢lnd⁡NN.

Istnieje bardzo dużo efektywnych konstrukcji ciągów punktów o niskiej dyskrepancji. Teraz poznamy jedynie kilka najbardziej popularnych z nich.

15.4.1. Ciąg Van der Corputa

Niech b≥2 będzie ustaloną liczbą naturalną. Wtedy dowolną liczbę naturalną n można jednoznacznie przedstawić w postaci

n⁢=⁢∑j=0∞aj⁢n⁢⁢bj,

gdzie aj∈0,1,…,b-1 i tylko dla skończonej liczby indeksów j mamy aj≠0. Mówiąc inaczej, aj są kolejnymi cyframi rozwinięcia liczby n w bazie b. Wykorzystując powyższe rozwinięcie funkcję radykalnej odwrotności ψb:0,1,2,…→0,1 definujemy jako

ψb⁢n⁢:=⁢∑j=1∞aj⁢n⁢⁢b-j+1.

Na przykład, dla bazy b=2 mamy 13=20+22+23=11012, czyli ψ2⁢13=12+18+116=1116.

Kolejne wartości radykalnej odwrotności,

ψb⁢1,ψb⁢2,…,ψb⁢n,…,

tworzą jednowymiarowy ciąg Van der Corputa. Dla b=2 kolejne punkty ciągu wynoszą więc

12,⁢14,⁢34,⁢18,⁢58,⁢38,⁢78,⁢116,⁢916,⁢516,⁢1316,⁢316,⁢1116,⁢716,⁢1516,⁢132,⁢1732,⁢…

Ciąg Van der Corputa jest ciągiem o niskiej dyskrepancji dla dowolnie dobranej bazy b, tzn. dyskrepancja N początkowych wyrazów szacuje się z góry przez C1⁢ln⁡N/N. Zauważmy jednak, że dla N=bk-1, k≥1, dostajemy równomierne rozmieszczenie punktów, tzn. tworzą one zbiór j/N+1: 1≤j≤N, którego dyskrepancja wynosi N+1-1. Stąd, pożądane są raczej mniejsze bazy b; im większe b tym rzadziej ze względu na N osiągana jest dyskrepancja proporcjonalna do 1/N.

15.4.2. Konstrukcje Haltona i Sobol'a

Jednowymiarowy ciąg Van der Corputa jest podstawą konstrukcji wielu ciągów w wyższych wymiarach d. Jedna z nich prowadzi do ciągu Haltona h→kk≥1, którego k-ty punkt wynosi

h→k⁢=⁢ψb1⁢k,ψb2⁢k,…,ψbd⁢kT.

Liczby b1,…,bd są tu danymi bazami. Od razu zauważamy, że wybór b1=⋯=bd nie jest dobry, bo prowadzi do siatki równomiernej w 0,1d, zob. rozdział 15.2. Jeśli jednak bi są liczbami pierwszymi to h→kk≥1 jest już ciągiem o niskiej dyskrepancji.

Minusem konstrukcji Haltona jest to, że dla dużych wymiarów d bazy bd też są duże co, jak wcześniej zauważyliśmy, nie jest korzystne. Problemu tego unika bardziej skomplikowana konstrukcja Sobol'a, gdzie na każdej zmiennej pracujemy z tą samą bazą b=2.

Ideę konstrukcji ciągu Sobol'a s→kk≥1 przedstawimy najpierw zakładając d=1. Niech

a→⁢k=a0⁢k,…,ar-1⁢kT

będzie wektorem kolejnych bitów w rozwinięciu dwójkowym liczby k,

k⁢=⁢a0⁢k+2⁢a1⁢k+⋯+2r-1⁢ar-1⁢k.

Wtedy

sk⁢=⁢y1⁢k2+y2⁢k4+⋯+yr⁢k2r,

gdzie

y1⁢ky2⁢k⋮yr⁢k⁢=⁢V⋅a0⁢ka1⁢k⋮ar-1⁢k⁢⁢mod⁢ 2,

a V jest specjalnie dobraną macierzą o elementach 0 i 1, zwaną macierzą generującą. (Zauważmy, że jeśli V jest identycznością to otrzymany ciąg jest ciągiem Van der Corputa.)

Oznaczając V=v→1,…,v→r możemy równoważnie zapisać

y→⁢k⁢=⁢a0⁢k⋅v→1⊕⋯⊕ar-1⁢k⋅v→r,

gdzie ⊕ jest operacją XOR, czyli dodawaniem bitów modulo 2,

0⊕0=0,⁢0⊕1=1,⁢1⊕0=1,⁢1⊕1=0.

Dla d≥2, kolejne współrzędne punktu s→k wyliczane są według powyższej recepty, ale używając różnych macierzy generujących V. I właśnie problem wyboru macierzy generujących tak, aby otrzymać ciąg o niskiej dyskrepancji jest istotą konstrukcji Sobol'a. Dodajmy, że jest to problem wysoce nietrywialny.

15.4.3. Sieci t,m,d i ciągi t,d

Dla b≥2 definiujemy b-prostokąt w 0,1d jako

a1bj1,a1-1bj1×⋯×adbjd,ad-1bjd,

gdzie ji∈0,1,2,…, ai∈0,1,…,bji-1, 1≤i≤d. Na przykład, jeśli d=1 i b=2 to przedziały 3/4,1 i 3/4,7/8 są b-prostokątami, ale nie jest nim 5/8,7/8. Zauważmy, że objętość b-prostokąta wynosi b-j1+⋯+jd.

Definicja 15.5

Niech b≥2 i 0≤t≤m. Siecią t,m,d w bazie b nazywamy zbiór bm punktów w 0,1d o własności, że w każdym b-prostokącie o objętości bt-m znajduje się dokładnie bt punktów tego zbioru.

Mówiąc inaczej, sieć t,m,d dokładnie pokazuje objętości b-prostokątów poprzez iloraz bt/bm liczby punktów, które należą do prostokąta do liczby wszystkich punktów.

Definicja 15.6

Ciąg nieskończony t→1,t→2,… w 0,1d nazywamy ciągiem t,d w bazie b jeśli dla wszystkich m>t i j=0,1,2,… segment

⁢t→i⁢:⁢j⁢bm<i≤j+1⁢bm⁢

jest siecią t,m,d w bazie b.

Następujące twierdzenie jest podstawą wielu konstrukcji ciągów o niskiej dyskrepancji.

Twierdzenie 15.2

Każdy ciąg t,d w bazie b jest ciągiem o niskiej dyskrepancji.

Pokazanie konkretnych konstrukcji ciągów t,d wykracza poza ramy tego wykładu. Powiemy tylko, że szczególne wybory macierzy generujących w konstrukcji Sobol'a prowadzą do ciągów t,d. Inne konstrukcje należą do Faure'a, Niederreitera, Tezuki i in.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.