Zagadnienia

15.1 Co to są metody quasi-Monte Carlo?
15.2 Dyskrepancja
15.3 Błąd quasi-Monte Carlo
- 15.3.1 Formuła Zaremby
- 15.3.2 Nierówność Koksmy-Hlawki
15.4 Ciągi o niskiej dyskrepancji

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)$ są $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.

	$\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}.$

Matematyka obliczeniowa II

Zagadnienia

15. Metody quasi-Monte Carlo

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

Definicja 15.1

15.2. Dyskrepancja

Definicja 15.2

15.3. Błąd quasi-Monte Carlo

15.3.1. Formuła Zaremby

Lemat 15.1

Dowód

15.3.2. Nierówność Koksmy-Hlawki

Definicja 15.3

Twierdzenie 15.1

Dowód

15.4. Ciągi o niskiej dyskrepancji

Definicja 15.4

15.4.1. Ciąg Van der Corputa

15.4.2. Konstrukcje Haltona i Sobol'a

15.4.3. Sieci i ciągi

Definicja 15.5

Definicja 15.6

Twierdzenie 15.2

15.4.3. Sieci $(t,m,d)$ i ciągi $(t,d)$