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 {\mathcal{F}}_{r}(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,

I_{d}(f)\,:=\,\int _{D}f(\vec{x})\, d\vec{x},\qquad D=[0,1]^{d}.
Definicja 15.1

Metoda quasi-Monte Carlo polega na przybliżeniu I_{d}(f) średnią arytmetyczną,

{\rm QMC}_{{d,N}}(f)\,=\,{\rm QMC}_{{d,N}}(\vec{t}_{1},\ldots,\vec{t}_{N})\,:=\,\frac{1}{n}\cdot\sum _{{j=1}}^{N}f(\vec{t}_{j}),

gdzie \vec{t}_{1},\vec{t}_{2},\ldots,\vec{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 \vec{x}=[x_{1},\ldots,x_{d}]^{T}\in D niech

[\vec{0},\vec{x})\,:=\,[0,x_{1})\times\cdots\times[0,x_{d})

oznacza d wymiarowy prostokąt w D, ,,zakotwiczony” w zerze. Zakładamy przy tym, że [\vec{0},\vec{0})=\emptyset.

Definicja 15.2

Dyskrepancją (z gwiazdką) danego zbioru N punktów \vec{t}_{j}\in[0,1)^{d}, 1\le j\le N, nazywamy wielkość

{\rm DISC}_{d}^{*}(\vec{t}_{1},\ldots,\vec{t}_{N})\,=\,\sup _{{\vec{x}\in D}}\,\left|\,{\rm DISC}_{d}(\vec{x};\vec{t}_{1},\ldots,\vec{t}_{N})\,\right|,

gdzie

{\rm DISC}_{d}(\vec{x};\vec{t}_{1},\ldots,\vec{t}_{N})\,=\, x_{1}\ldots x_{d}-\frac{\#\left\{ j\,:\;\vec{t}_{j}\in[\vec{0},\vec{x})\right\}}{N},

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

Ponieważ x_{1}\cdots x_{d} jest d-wymiarową objętością prostokąta [\vec{0},\vec{x}), dyskrepancja pokazuje jak dobrze danych N punktów aproksymuje objętości tych prostokątów. Równoważnie, oznaczając przez {\bf 1}_{{[\vec{0},\vec{x})}} funkcję charakterystyczną prostokąta [\vec{0},\vec{x}) mamy

\int _{D}{\bf 1}_{{[\vec{0},\vec{x})}}(t)\, dt=x_{1}\cdots x_{d}\qquad\mbox{oraz}\qquad\frac{1}{N}\sum _{{j=1}}^{N}{\bf 1}_{{[\vec{0},\vec{x})}}(\vec{t}_{j})=\frac{\#\{ j\,:\;\vec{t}_{j}\in[\vec{0},\vec{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 t_{j}

{\rm DISC}_{1}^{*}(t_{1},\ldots,t_{N})\,\ge\,\frac{1}{2N}. (15.1)

Rzeczywiście, {\rm DISC}_{1}(x;t_{1},\ldots,t_{N}) jako funkcja x ma na każdym z przedziałów [0,t_{1}), [t_{{j-1}},t_{j}), 2\le j\le N, pochodną równą 1, oraz przyjmuje zero dla x=0. Zatem, jeśli dyskrepancja byłaby mniejsza od 1/(2N) to mielibyśmy

t_{1}<\frac{1}{2N},\qquad t_{j}-t_{{j-1}}<\frac{1}{N},\quad 2\le j\le N,

a stąd

t_{N}\,=\, t_{1}+\sum _{{j=2}}^{N}(t_{j}-t_{{j-1}})\,<\,\frac{1}{2N}+\frac{N-1}{N}\,=\, 1-\frac{1}{2N}.

Otrzymujemy sprzeczność, bo dla t_{N}<x<1-1/(2N)

{\rm DISC}_{1}(x;t_{1},\ldots,t_{N})\,<\,\left(1-\frac{1}{2N}\right)-1\,=\,-\frac{1}{2N}.

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

t_{j}^{*}\,=\,\frac{j-1/2}{N},\qquad 1\le j\le N.

W tym przypadku algorytm {\rm QMC}_{{1,N}} redukuje się do zasady punktu środkowego.

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

t_{1},t_{2},\ldots,t_{n},\ldots\,\subset[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ść

{\rm DISC}_{1}^{*}(t_{1},\ldots,t_{N})\,\ge\, c\,\frac{\ln N}{N}

zachodzi dla nieskończenie wielu N.

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

\vec{t}_{{j_{1},\ldots,j_{d}}}\,=\,[t_{{j_{1}}}^{*},\ldots,t_{{j_{d}}}^{*}]^{T},\qquad 1\le j_{i}\le n,1\le i\le d,

będzie równomiernym rozkładem N=n^{d} punktów w D. Rozpatrzenie prostokąta

[0,\tfrac{1}{2n})\times\underbrace{[0,1)\times\cdots\times[0,1)}_{{d-1}}

wystarcza, aby się przekonać, że wtedy dyskrepancja wynosi co najmniej \tfrac{1}{2n}=\tfrac{1}{2}N^{{-1/d}}.

15.3. Błąd quasi-Monte Carlo

15.3.1. Formuła Zaremby

Jeśli f:[0,1]\to\mathbb{R} jest funkcją różniczkowalną to dla każdego x mamy

f(x)\,=\, f(1)-\int _{x}^{1}f^{{\prime}}(t)\, dt\,=\, f(1)-\int _{0}^{1}{\bf 1}_{{(x,1]}}(t)\, f^{{\prime}}(t)\, dt. (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,

\emptyset\,\ne\, U\,\subseteq\,\{ 1,2,\ldots,d\},

niech D_{U}=[0,1]^{{|U|}}, gdzie |U| jest liczbą elementów w U. Niech dalej \vec{x}_{U}\in D_{U} będzie wektorem powstałym z wektora \vec{x}=[x_{1},x_{2},\ldots,x_{d}]^{T}\in D poprzez usunięcie współrzędnych x_{j} z j\notin U, a (\vec{x}_{U};1)\in D wektorem, którego j-ta współrzędna wynosi x_{j} dla j\in U oraz wynosi 1 dla j\notin U. Na przykład, jeśli d=5 i U=\{ 1,4\} to dla \vec{x}=[x_{1},x_{2},x_{3},x_{4},x_{5}]^{T} mamy \vec{x}_{U}=[x_{1},x_{4}]^{T} i (\vec{x}_{U};1)=[x_{1},1,1,x_{4},1]^{T}.

W końcu, niech

f_{U}^{{\prime}}\,=\,\frac{\partial^{{|U|}}f}{\prod _{{j\in U}}\partial u_{j}}

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

Lemat 15.1

Jeśli funkcja f:D\to\mathbb{R} ma ciągłe pochodne cząstkowe mieszane f_{U}^{{\prime}} dla wszystkich \emptyset\ne U\subseteq\{ 1,2,\ldots,d\} to

f(\vec{x})\,=\, f(\vec{1})\,+\,\sum _{U}(-1)^{{|U|}}\int _{{D_{U}}}{\bf 1}_{{(\vec{x}_{U},\vec{1}_{U}]}}(\vec{z}_{U})\, f_{U}^{{\prime}}(\vec{z}_{U};1)\, d\vec{z}_{U},\qquad\vec{x}\in D (15.3)

(\vec{1}=[1,1,\ldots,1]\in\mathbb{R}^{d}).

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\ge 2. Stosując założenie indukcyjne do f z ustaloną ostatnią współrzędną x_{d} mamy

f(\vec{x})\;=\; f(\vec{x}_{{\{ d\}}};1)\,+\,\sum _{V}(-1)^{{|V|}}\int _{{D_{V}}}{\bf 1}_{{(\vec{x}_{V},\vec{1}_{V}]}}(\vec{z}_{V})\, f_{V}^{{\prime}}(\vec{z}_{V},x_{d};1)\, d\vec{z}_{V},

gdzie sumowanie jest po wszystkich \emptyset\ne V\subseteq\{ 1,2,\ldots,d-1\}. Stosując dalej wzór (15.2) ze względu na x_{d} odpowiednio do f(\vec{x}_{{\{ d\}}};1) i f_{V}^{{\prime}}(\vec{z}_{V},x_{d};1) dostajemy

\displaystyle f(\vec{x}) \displaystyle= \displaystyle f(\vec{1})\,-\,\int _{{D_{{\{ d\}}}}}{\bf 1}_{{\left(\vec{x}_{{\{ d\}}},\vec{1}_{{\{ d\}}}\right]}}\left(\vec{z}_{{\{ d\}}}\right)\, f_{{\{ d\}}}^{{\prime}}\left(\vec{z}_{{\{ d\}}};1\right)\, d\vec{z}_{{\{ d\}}}
\displaystyle+(-1)^{{|V|}}\int _{{D_{V}}}{\bf 1}_{{(\vec{x}_{V},\vec{1}_{V}]}}(\vec{z}_{V})\, f_{V}^{{\prime}}(\vec{z}_{V};1)\, d\vec{z}_{V}
\displaystyle+\sum _{{V\cup\{ d\}}}(-1)^{{|V\cup\{ d\}|}}\int _{{D_{{V\cup\{ d\}}}}}{\bf 1}_{{\left(\vec{x}_{{V\cup\{ d\}}},\vec{1}_{{V\cup\{ d\}}}\right]}}\left(\vec{z}_{{V\cup\{ d\}}}\right)\, f_{{V\cup\{ d\}}}^{{\prime}}\left(\vec{z}_{{V\cup\{ d\}}};1\right)\, d\vec{z}_{{V\cup\{ d\}}}.

Teraz wystarczy porównać otrzymaną formułę do prawej strony (15.3). Drugi składnik odpowiada U=\{ d\}, trzeci składnik podzbiorom U\ne\emptyset do których nie należy d, a czwarty podzbiorom U takim, że d\in U i |U|\ge 2.

Zauważmy, że rozwijając f(\vec{x}) i f(\vec{t}_{j}) zgodnie ze wzorem (15.3) mamy

\int _{D}f(\vec{x})\, d\vec{x}\,=\, f(\vec{1})+\sum _{U}(-1)^{{|U|}}\int _{{D_{U}}}\left(\int _{D}{\bf 1}_{{(\vec{x}_{U},\vec{1}_{U}]}}(\vec{z}_{U})\, d\vec{x}\right)\, f_{U}^{{\prime}}(\vec{z}_{U};1)\, d\vec{z}_{U},
\frac{1}{N}\,\sum _{{j=1}}^{N}f(\vec{t}_{j})\,=\, f(\vec{1})+\sum _{U}(-1)^{{|U|}}\int _{{D_{U}}}\left(\frac{1}{N}\sum _{{j=1}}^{N}{\bf 1}_{{((\vec{t}_{j})_{U},\vec{1}_{U}]}}(\vec{z}_{U})\right)\, f_{U}^{{\prime}}(\vec{z}_{U};1)\, d\vec{z}_{U}.

Ponieważ wartość funkcji charakterystycznej odcinka (\vec{a},\vec{b}] w punkcie \vec{c} jest równa wartości funkcji charakterystycznej odcinka [\vec{0},\vec{c}) w punkcie \vec{a} to

\int _{D}{\bf 1}_{{(\vec{x}_{U},\vec{1}_{U}]}}(\vec{z}_{U})\, d\vec{x}\,=\,\int _{{D_{U}}}{\bf 1}_{{[\vec{0}_{U},\vec{z}_{U})}}(\vec{x}_{U})\, d\vec{x}_{U}\,=\,\prod _{{j\in U}}z_{j},
\frac{1}{N}\,\sum _{{j=1}}^{N}{\bf 1}_{{((\vec{t}_{j})_{U},\vec{1}_{U}]}}(\vec{z}_{U})\,=\,\frac{1}{N}\,\#\left\{ j\,:\;(\vec{t}_{j})_{U}\in[\vec{0}_{U},\vec{z}_{U})\right\}\,=\,\frac{1}{N}\,\#\left\{ j\,:\;\vec{t}_{j}\in\left[\vec{0},(\vec{z}_{U};1)\right)\right\}.

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

I_{d}(f)-{\rm QMC}_{{d,N}}(f)\,=\,\sum _{U}(-1)^{{|U|}}\int _{{D_{U}}}{\rm DISC}_{d}\left((\vec{z}_{U};1);\vec{t}_{1},\ldots,\vec{t}_{N}\right)\cdot f_{U}^{{\prime}}(\vec{z}_{U};1)\, d\vec{z}_{U}. (15.4)

15.3.2. Nierówność Koksmy-Hlawki

Oznaczmy przez {\mathcal{V}}_{d}(D) klasę funkcji f:D\to\mathbb{R}, których pochodne mieszane f_{U}^{{\prime}} istnieją i są ciągłe dla wszystkich \emptyset\ne U\subseteq\{ 1,\ldots,d\}.

Definicja 15.3

Wahaniem (w sensie Hardy-Krause) funkcji f\in{\mathcal{V}}_{d}(D) nazywamy wielkość

V_{d}(f)\,:=\,\sum _{U}\int _{{D_{U}}}\left|f_{U}^{{\prime}}(\vec{z}_{U};1)\right|\, d\vec{z}_{U}.

Zauważmy, że dla d=1,

V_{1}(f)\,=\,\int _{0}^{1}|f^{{\prime}}(x)|\, dx\,=\,\sup\left\{\,\sum _{{j=1}}^{k}|f(z_{j})-f(z_{{j-1}})|\,:\; k\ge 1,0=z_{0}<z_{1}<\ldots<z_{k}=1\,\right\}

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\in{\mathcal{V}}_{d}(D) to błąd metody quasi-Monte Carlo

{\rm QMC}_{{d,N}}(f)=\frac{1}{N}\sum _{{j=1}}^{N}f(\vec{t}_{j})

dla aproksymacji całki I_{d}=\int _{D}f(\vec{x})\, d\vec{x} szacuje się przez

\left|I_{d}(f)-{\rm QMC}_{{d,N}}(f)\right|\;\le\;{\rm DISC}^{*}_{d}(\vec{t}_{1},\ldots,\vec{t}_{N})\cdot V_{d}(f).
Dowód

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

\displaystyle|I_{d}(f)-{\rm QMC}_{{d,N}}(f)| \displaystyle\le \displaystyle\sum _{U}\int _{{D_{U}}}\left|{\rm DISC}_{d}\left((\vec{z}_{U};1);\vec{t}_{1},\ldots,\vec{t}_{N}\right)\right|\cdot\left|f_{U}^{{\prime}}(\vec{z}_{U};1)\right|\, d\vec{z}_{U}
\displaystyle\le \displaystyle{\rm DISC}_{d}^{*}(\vec{t}_{1},\ldots,\vec{t}_{N})\cdot\sum _{U}\int _{{D_{U}}}\left|f_{U}^{{\prime}}(\vec{z}_{U};1)\right|\, d\vec{z}_{U}.

W nierówności Koksmy-Hlawki czynnik błędu V_{d}(f) zależny jedynie od funkcji jest oddzielony od czynnika błędu {\rm DISC}_{d}^{*}(\vec{t}_{1},\ldots,\vec{t}_{N}) zależnego od wyboru punktów. O ile nie mamy wpływu na wahanie funkcji, możemy starać się wybrać punkty \vec{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 {\mathcal{V}}_{d}(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\in{\mathcal{V}}_{d}(D) z V(f)\le 1 jedynie do oszacowania z dołu błędu przez cN^{{-1}} (a nie cN^{{-r/d}}).

Okazuje się, że pełne odpowiedzi na zadane pytania nie są znane. Najlepsze ciągi nieskończone \vec{t}^{{\,\,}}_{1},\vec{t}^{{\,\,}}_{2},\ldots,\vec{t}^{{\,\,}}_{n},\ldots spełniają nierówność

{\rm DISC}_{d}^{*}\left(\vec{t}_{1}^{{\,\,}},\ldots,\vec{t}_{N}^{{\,\,}}\right)\,\le\, C_{d}\,\frac{\ln^{d}N}{N},\qquad N=1,2,3,\ldots,

gdzie C_{d}>0 nie zależy od N.

Wydaje się więc, że w klasie {\mathcal{V}}_{d}(D) metody quasi-Monte Carlo dają istotnie lepsze wyniki od Monte Carlo, bo błąd nie tylko jest deterministyczny, ale też dla f\in{\mathcal{V}}_{d}(D) zbiega do zera dużo szybciej, tzn. błąd jest rzędu N^{{-1+\epsilon}} dla dowolnego \epsilon>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 C_{d}\ln^{d}N może mieć istotne znaczenie i wręcz sprawiać, że nierówność Koksmy-Hlawki staje się bezużyteczna. Dodatkowo, dobre oszacowanie C_{d} 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 \vec{t}_{1},\vec{t}_{2},\vec{t}_{3},\ldots nazywamy ciągem o niskiej dyskrepancji jeśli istnieje C_{d}>0 taka, że dla wszystkich N

{\rm DISC}_{d}^{*}\left(\vec{t}_{1},\ldots,\vec{t}_{N}\right)\,\le\, C_{d}\,\frac{\ln^{d}N}{N}.

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\ge 2 będzie ustaloną liczbą naturalną. Wtedy dowolną liczbę naturalną n można jednoznacznie przedstawić w postaci

n\,=\,\sum _{{j=0}}^{\infty}a_{j}(n)\, b^{j},

gdzie a_{j}\in\{ 0,1,\ldots,b-1\} i tylko dla skończonej liczby indeksów j mamy a_{j}\ne 0. Mówiąc inaczej, a_{j} są kolejnymi cyframi rozwinięcia liczby n w bazie b. Wykorzystując powyższe rozwinięcie funkcję radykalnej odwrotności \psi _{b}:\{ 0,1,2,\ldots\}\to[0,1) definujemy jako

\psi _{b}(n)\,:=\,\sum _{{j=1}}^{\infty}a_{j}(n)\, b^{{-(j+1)}}.

Na przykład, dla bazy b=2 mamy 13=2^{0}+2^{2}+2^{3}=(1101)_{2}, czyli \psi _{2}(13)=\tfrac{1}{2}+\tfrac{1}{8}+\tfrac{1}{16}=\tfrac{11}{16}.

Kolejne wartości radykalnej odwrotności,

\psi _{b}(1),\psi _{b}(2),\ldots,\psi _{b}(n),\ldots,

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

\frac{1}{2},\,\frac{1}{4},\,\frac{3}{4},\,\frac{1}{8},\,\frac{5}{8},\,\frac{3}{8},\,\frac{7}{8},\,\frac{1}{16},\,\frac{9}{16},\,\frac{5}{16},\,\frac{13}{16},\,\frac{3}{16},\,\frac{11}{16},\,\frac{7}{16},\,\frac{15}{16},\,\frac{1}{32},\,\frac{17}{32},\,\ldots

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 C_{1}\ln N/N. Zauważmy jednak, że dla N=b^{k}-1, k\ge 1, dostajemy równomierne rozmieszczenie punktów, tzn. tworzą one zbiór \{ j/(N+1):\, 1\le j\le 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 \{\vec{h}_{k}\} _{{k\ge 1}}, którego k-ty punkt wynosi

\vec{h}_{k}\,=\,\left[\psi _{{b_{1}}}(k),\psi _{{b_{2}}}(k),\ldots,\psi _{{b_{d}}}(k)\right]^{T}.

Liczby b_{1},\ldots,b_{d} są tu danymi bazami. Od razu zauważamy, że wybór b_{1}=\cdots=b_{d} nie jest dobry, bo prowadzi do siatki równomiernej w [0,1)^{d}, zob. rozdział 15.2. Jeśli jednak b_{i} są liczbami pierwszymi to \{\vec{h}_{k}\} _{{k\ge 1}} jest już ciągiem o niskiej dyskrepancji.

Minusem konstrukcji Haltona jest to, że dla dużych wymiarów d bazy b_{d} 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 \{\vec{s}_{k}\} _{{k\ge 1}} przedstawimy najpierw zakładając d=1. Niech

\vec{a}(k)=[a_{0}(k),\ldots,a_{{r-1}}(k)]^{T}

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

k\,=\, a_{0}(k)+2a_{1}(k)+\cdots+2^{{r-1}}a_{{r-1}}(k).

Wtedy

s_{k}\,=\,\frac{y_{1}(k)}{2}+\frac{y_{2}(k)}{4}+\cdots+\frac{y_{r}(k)}{2^{r}},

gdzie

\left[\begin{array}[]{c}y_{1}(k)\\
y_{2}(k)\\
\vdots\\
y_{r}(k)\end{array}\right]\,=\, V\cdot\left[\begin{array}[]{c}a_{0}(k)\\
a_{1}(k)\\
\vdots\\
a_{{r-1}}(k)\end{array}\right]\;\mbox{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=[\vec{v}_{1},\ldots,\vec{v}_{r}] możemy równoważnie zapisać

\vec{y}(k)\,=\, a_{0}(k)\cdot\vec{v}_{1}\oplus\cdots\oplus a_{{r-1}}(k)\cdot\vec{v}_{r},

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

0\oplus 0=0,\quad 0\oplus 1=1,\quad 1\oplus 0=1,\quad 1\oplus 1=0.

Dla d\ge 2, kolejne współrzędne punktu \vec{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\ge 2 definiujemy b-prostokąt w [0,1)^{d} jako

\left[\frac{a_{1}}{b^{{j_{1}}}},\frac{a_{1}-1}{b^{{j_{1}}}}\right)\times\cdots\times\left[\frac{a_{d}}{b^{{j_{d}}}},\frac{a_{d}-1}{b^{{j_{d}}}}\right),

gdzie j_{i}\in\{ 0,1,2,\ldots\}, a_{i}\in\{ 0,1,\ldots,b^{{j_{i}-1}}\}, 1\le i\le d. Na przykład, jeśli d=1 i b=2 to przedziały [3/4,1) i [3/4,7/8)b-prostokątami, ale nie jest nim [5/8,7/8). Zauważmy, że objętość b-prostokąta wynosi b^{{-(j_{1}+\cdots+j_{d})}}.

Definicja 15.5

Niech b\ge 2 i 0\le t\le m. Siecią (t,m,d) w bazie b nazywamy zbiór b^{m} punktów w [0,1)^{d} o własności, że w każdym b-prostokącie o objętości b^{{t-m}} znajduje się dokładnie b^{t} punktów tego zbioru.

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

Definicja 15.6

Ciąg nieskończony \vec{t}_{1},\vec{t}_{2},\ldots w [0,1)^{d} nazywamy ciągiem (t,d) w bazie b jeśli dla wszystkich m>t i j=0,1,2,\ldots segment

\left\{\,\vec{t}_{i}\,:\; jb^{m}<i\le(j+1)b^{m}\,\right\}

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.