Seminar Numerik stochastischer Modelle

Low-Discrepancy Folgen (Low-Discrepancy Sequences)und die quasi Monte Carlo Methode

Mathematisch-Naturwissenschaftlichen Fakultät IHumboldt-Universität zu Berlin

vonHermann Schwarz

Inhaltsverzeichnis

1 Einführung 1

2 Koksma-Hlawka Ungleichung [3] 2

3 Folgen kleiner Diskrepanz 63.1 Halton Folge [3] . . . . . . . . . . . . . . . . . . . . . . . . . . 63.2 Van der Corput, Basis 2 . . . . . . . . . . . . . . . . . . . . . 73.3 Faure Folge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73.4 Simulation [3] . . . . . . . . . . . . . . . . . . . . . . . . . . . 93.5 Vergleich der Monte Carlo und Quasi-Monte Carlo Methoden . 9

i

Kapitel 1

Einführung

Als erstes möchte ich, dass wir uns an die Monte Carlo Mathode erinnern.Es handelt sich hierbei um eine Methode zur numerischen Integration mitHilfe einer Folge von Zufallszahlen:

1∫0

f(u)du ≈ 1

N

N∑i=1

f(ξi) (1.1)

Die Funktion f ist eine bekannte Funktion, deren Integral wir approximie-ren möchten. Dafür wird die Summe der Funktionswerte dieser Funktion inden zufällig gewählten Punkten, multipliziert mit einem N-tel, berechnet. Nist dabei die Anzahl der zufälligen Punkte, und ξ1, ..., ξN sind diese zufälliggewählte Punkte aus dem Intervall [0, 1].Wenn die Punkte mit Hilfe folgender Regel gewählt werden: ξi = i

Nheißt

die Methode Trapezregel. Wenn die Punkte randomisiert oder pseudorando-misiert ausgewählt werden, handelt es sich um die Monte Carlo Methode.Wenn die Punkte aber aus der sogenannten Low-Discrepancy Folgen ausge-wählt werden, wird die Methode als quasi-Monte Carlo Methode bezeichnet.Es wird heute also wieder um die Monte Carlo Methode gehen, allerdings ineiner speziellen Form.Welche Motivation steckt hinter der quasi-Monte Carlo Methode? Die Folgenkleiner Diskrepanz, die auch quasi-random Folgen genannt werden, erlaubenin vielen Fällen eine Verbesserung der Performanz der Monte Carlo Methode,also Kürzere Rechenzeiten und auch größere Genauigkeit. Die Genauigkeitder Berechnung des Integrals wird dadurch erhöht, dass eine Folge kleinerDiskrepanz gleichmäßiger über das Intervall zerstreut ist, über dem Integralberechnet werden soll.

1

Kapitel 2

Koksma-Hlawka Ungleichung [3]

Diese Ungleichung ist ein bemerkenswertes Ergebnis, das uns zeigt, dass derFehler, der bei der Verwendung der quasi Monte Carlo Methoden auftrittbegrenzt ist, und zwar durch das Produkt zweier Terme. Einer der beidenTerme ist von der Funktion f abhängig, deren Intergral wir mit Hilfe der qua-si Monte Carlo Methode approximieren. Der andere Term ist die Diskrepanzder Menge ξ1, ..., ξN .Bevor ich die Ungelichung zeige, muss definiert werden, was eine Funktionendlicher Variation ist.

NotationFür x ∈ [0, 1]s bezeichnen wir mit x̌ den Punkt, der sich symmetrisch bzgl.zum Zentrum vom Hyperwürfel [0, 1]s. Wenn eine Funktion auf diesem Hyper-würfel definiert ist, f : [0, 1]s → R, bezeichnen wir mit f̌ den Funktionswertdieser Funktion im Punkt x̌ (f̌(x) = f(x̌)).

Definition (Funktionen endlicher Variation)Eine Funktion f : [0, 1]s → R besitzt endliche Variation, falls ein Maß µauf der σ-Algebra der Borelmengen des Intervalls [0, 1]s mit dem Träger in[0, 1]s \ 0 existiert, so dassf̌(x) = f̌(0) + µ[[0, x]] ∀x ∈ [0, 1]s

Die Norm dieses Masses ||µ|| wird als Variation von f bezeichnet und wirdmit V (f). Da ich den Begriff der Diskrepanz einer Folge mehrmals verwendenwerde, möchte ich ihn nochmal formulieren:

2

3

Definition (Diskrepanz) [2]Diskrepanz ist ein Maß dafür, wie weit die durch Folge ξ1, ..., ξN gegebeneempirische Verteilung von der Gleichverteilung auf [0, 1] entfernt ist.Formal:

D∗∞(ξ, N) = sup

R∈J∗|A(R; ξ)

N− λs(R)| (2.1)

wo ξ die Folge ξ1, ..., ξN ist, λs ist s-dimensionales Lebesgue Maß, A(R; ξ) istdie Anzahl der Punkte aus ξ, die in R liegen.J∗ ist die Menge der Intervalle der Form:s∏

i=1

[0, ui) , wobei ui im halboffenen Intervall [0, 1) liegt.

Definition (Koksma-Hlawka Ungleichung) [3]Sei f eine Funktion mit endlicher Variation auf [0, 1]s und sei ξ = (ξn) eineFolge in [0, 1]s, dann gilt für irgendein N∣∣∣∣∣∣∣

∫[0,1]s

f(u)du− 1

N

N∑n=1

f(ξn)

∣∣∣∣∣∣∣ ≤ V (f) ·D∗∞(ξ, N) (2.2)

Bevor wir diesen Satz Beweisen, müssen wir folgendes Lemma voraussetzen:Wir Bezeichnen mit µ̌ die Abbildung µ(x̌), wobei x̌ der bezüglich des Zen-trums des Einheitswürfels [0, 1]s symmetrische Punkt ist.LemmaSei µ ein Maß auf [0, 1]s, mit Verteilungsfunktion F . Dann gilt mit der Be-zeichnung aus dem ersten Vortrag:∫

π(x)dµ̌(x) =

∫F (x)dx (2.3)

Der Beweis der Koksma-Hlawka Ungleichung:Wir bezeichnen F := f̌ − f̌(0)Dann

1

N

N∑n=1

f(ξn)−∫

f(x)dx =1

N

N∑n=1

f̌(ξ̌n)−∫

f(x̌)dx =

1

N

N∑n=1

F (ξ̌n)−∫

F (x̌)dx =

∫ (1

N

N∑n=1

1x̌≤ξ̌n− π(x)

)dµ̌(x) =

4

∫ (1

N

N∑n=1

1ξn≤x − π(x)

)dµ̌(x)

Daraus folgt:∣∣∣∣∣ 1

N

N∑n=1

f(ξn)−∫

f(x)dx

∣∣∣∣∣ ≤ D∗∞(ξ, N)||µ̌|| = D∗

∞(ξ, N)||µ||

Diese Ungleichung stellt die Basis für die Überlegenheit der Quasi-MonteCarlo gegenüber der Monte Carlo Methode, weil der Fehler der Berechnungdes Integrals gegen Null strebt, und zwar mit Geschwindigkeit O( (logN)s

N).

Allerdings, wie in einem Paper von zwei Japanischen Mathematikern HozumiMorohoshi und Masanori Fushimi über Fehlerabschätzung von quasi-MonteCarlo Methode beschrieben wurde, ist es fast unmöglich, Fehlerabschätzungmit Hilfe dieser Ungleichung zu machen, weil die Berechnung der Variationder Funktion f normalerweise unmöglich ist. Für die Fehlerabschätzung derquasi-Monte Carlo Methode gibt es deswegen Methoden, die für die Praxismehr geeignet sind [1].)

Für den Fall, wenn es doch möglich ist, die Variation der Funktion auszu-rechnen, kann man zur Erleichterung der Berechnung folgende Ungleichungenbenutzen [3]. Sie können versuchen, diese Ungleichungen zu beweisen:Seien f, g : [0, 1]s → R, und sei λ ∈ RDann gilt:

V (f + g) ≤ V (f) + V (g)V (λf) = |λ|V (f)||f ||∞ ≤ |f(1)|+ V (f)

Man kann auch folgende Eigenschaft der Variation benutzen:Seien f : [0, 1]s → R, f ∈ Cs

Dann gilt:

V (f) =

∣∣∣∣∣∣∣∣ δsf

δx1....δxs

∣∣∣∣∣∣∣∣1

An der Ungleichung sieht man, dass die Größe des Fehlers von der Diskre-panz abhängt. Um den Fehler zu minimieren ist es also sinnvoll, Mengen vonPunkten zu konstruieren, die möglichst gleichverteilt sind, d.h. mit möglichstkleiner Diskrepanz. Allerdings ist es unmöglich, die Diskrepanz beliebig klein

5

zu halten. Sie hat eine untere Schranke, die vom einem britischen Mathema-tiker Klaus Friedrich Roth entdeckt wurde:

Theorem von Roth [3]Sei s ≥ 2 und sei ξ eine Folge aus dem Intervall [0, 1]s.∀N ≥ 1 gilt:

D∗∞(ξ, N) ≥ D∗

2(ξ, N) > C(s)(logN)

(s−1)2

N(2.4)

wo C(s) > 0 eine Konstante ist, die nur von s abhängt.

Falls s ≥ 1∀N ≥ 1 gilt:

D∗∞(ξ, N) > C1(s)

(logN)(s)2

N(2.5)

BemerkungFür die Logaritmen zur Basis 2, die beiden Ungleichungen sind gültig mitC(s) = (22s+4(s− 1)

s−12 )−1 und

C1(s) = C(s + 1)− ε

Kapitel 3

Folgen kleiner Diskrepanz

3.1 Halton Folge [3]Sei p ≥ 2 eine natürliche Zahl. Wenn nunn = a0 + a1p + ... + akp

k die p-adische Entwicklung von n ∈ N ist,wobei ai ∈ {0, ...., p− 1}, i = 0, ...., k (k+1 ist die Anzahl der Stellen derbenutzten Zahl unter benutzten Basis). Dann ist die Zahl ϕp(n), eine Kom-ponente eines Elements der Halton Folge, folgendermassen definiert:ϕp(n) = a0

p+ a1

p2 + ... + ak

pk+1

Die Halton Folge x1, ..., xN mit xn ∈ [0, 1]s ist dann definiert durch:xn = (ϕp1(n), ....., ϕps(n))Die Zahlen p1 bis ps sind natürliche Zahlen größer Eins und paarweise relativprim.Wir erhalten also N Folgenglieder, die s-dimensionale Vektore sind.

Beispiel:wir möchten das 15. Glied der Halton Folge in einem 3-dimensionalen Raumausrechnen.Wählen wir als erstes die Zahlen p1, p2 und p3 aus. Da sie paarweise relativprim sein müssen, können wir z.B. 2,3 und 5 nehmen.p-adische Entwicklungen von 15 mit Basen 2,3 und 5:

Basis 2: 15 = 1 · 20 + 1 · 21 + 1 · 22 + 1 · 23 = 11112

Basis 3: 15 = 0 · 30 + 2 · 31 + 1 · 32 = 1203

Basis 5: 15 = 0 · 50 + 3 · 51 = 305

6

7

Nun berechnen wir die drei Komponenten des 15. Gliedes der Halton Fol-ge:ϕ2(15) = 1

2+ 1

22 + 123 + 1

24 = 0.11112 = 1516

ϕ3(15) = 03

+ 232 + 1

33 = 0.0213 = 727

ϕ5(15) = 05

+ 352 = 0.035 = 3

25

Diese drei Komponenten kann man auch als Spiegelung der Darstellung vonder Zahl n (in unserem Fall die 15) in Basis p (in unserem Beispiel 2,3 und5) am Komma bezeichnen.Das 15. Glied der Haltonfolge:x15 = (ϕp1(15), ϕp2(15), ϕp3(15)) = (ϕ2(15), ϕ3(15), ϕ5(15)) = (15

16, 7

27, 3

25)

Die Beschränkung der Diskrepanz:

D∗∞(x, N) ≤ 1

N

(s∑

i=1

pi + (2 log N)s

s∏i=1

(pi − 1)

logpi

)(3.1)

Kurz geschrieben: D∗∞(x, N) = O( (log N)s

N)

Was man vielleicht noch zu Halton Folge sagen könnte, dass bei den Fol-gen ab der Dimension 14 klare Strukturen zu erkennen sind.

3.2 Van der Corput, Basis 2Die eindimensionale Haltonfolge mit Basis 2 hat einen eigenen Namen, vander Corput Folge und gilt als die elementare Folge kleiner Diskrepanz undwird als Baustein für viele andere Folgen kleiner Diskrepanz verwendet. Dadie Folge eindimensional ist, möchte ich mit Hilfe dieser Folge veranschauli-chen, nach welchem Prinzip sie und die Halton Folge sich in einem Intervallbzw. in einem s-dimensionalen Würfel verteilen.(an der Tafel die ersten 8 Elemente der folge zeigen, als eine Übung anbieten,die ersten 8 Elemente auszurechnen )

3.3 Faure FolgeBei der Erzeugung der Faure Folgen, wie auch bei vielen anderen, wie z.B.Hammerslay Folgen, wird die Spiegelung der Darstellung von der Zahl n in

8

Basis p benutzt, wenn wir das n-te Glied einer Folge kleiner Diskrepanz aus-rechnen möchten.Faure Folgen sind also den Halton Folgen sehr ähnlich, und haben mit denennur zwei Unterschiede. Erstens für die Konstruktion dieser Folge wird nur ei-ne Basis für alle Dimensionen verwendet. Und zweitens werden hier für jedeDimension durch Permutationen gebildet.Sei s ≥ 1 und wir wählen eine Primzahl p ≥ s. Sei ∆p Menge der p-adischenEntwicklungen rationaler Zahlen aus dem Intervall [0, 1[Wir definieren eine Abbildung Cp : ∆p → ∆p, die die Menge der p-adischenEntwicklungen rationaler Zahlen auf sich selbst abbildet.Hier ist noch einmal die p-adische Entwicklung der rationalen Zahl n:ϕpi

(n) = a0

pi+ a1

p2i

+ ... + ak

pk+1i

Diese Entwicklung, wie wir schon wissen, ist die Darstellung der i-ten Kom-ponente eines Elements der Halton Folge.

Cp(ϕp(n)) =

kPi=0

(i0)ai mod p

p+

kPi=1

(i1)ai mod p

p2 + ... +

kPi=k

(ik)ai mod p

pk+1

Wenn wir die Abbildung Ci j-mal Anwenden, setzen wir sie gleich Null,wir schreiben Cj

i = 0, falls j > i ≥ 0Falls i ≥ j ≥ 0, Cj

i =(

ij

)Die Abbildung Cp ist injektiv und es gilt: (Cp)

p = Id, wenn man also dieAbbildung p-mal anwendet, bildet man die p-adische Entwicklung einer ra-tionaler Zahl auf sich selbst ab.Es gilt auch (Cp)

−1 = (Cp)p−1.

Sei nun s ≥ 2 und p ≥ 3 ist eine Primzahl, so dass p ≥ s.Die Faure Folge x1, ..., xN mit xn ∈ [0, 1]s ist dann definiert durch: xn =(ϕp(n), Cp(ϕp(n)), C2

p(ϕp(n)), ....., Cs−1p (ϕp(n)))

Die Basis der Faure Folge ist also die kleinste Primzahl, die größer gleich derDimensionzahl ist. Das erweist sich als vorteilhaft bei Problemen mit großerDimensionszahl. Z.B. wenn das Problem Dimensionsgröße 50 hat, würden wirfür die 50. Dimension eines Elements der Halton Folge die 50. Primzahl ver-wenden, nämlich 229. Für die Faure Folge müssten wir aber an dieser Stelledie Primzahl 53 nehmen. Damit werden die leeren Stellen mit Faure besserals mit Halton Folgen ausgefüllt, wenn das Problem hochdimensional ist.

Die Diskrepanz wird vom folgenden Term begrenzt:

D∗∞(x, N) ≤ 1

N

(1

s!

(p− 1

2 log p

)s

(log N)s + O((log N)s−1)

)(3.2)

9

3.4 Simulation [3]Angenommen, X ist eine gleichverteilte Zufallsvariable auf [0, 1]s und f isteine Riemann-integrierbare Funktion definiert auf disem s-dimensionalen In-tervall [0, 1]s. Anstatt nun den Erwartungswert folgendermassen zu berech-nen:

E [f(X)] = limN→∞

1

N

N∑n=1

f(Un)f.s. (3.3)

wo (Un) Folge von Borelmengen in [0, 1]s, können wir nun (Un) durch diegleichverteilte Folge (xn) in [0, 1]s ersetzen.Der Vorteil dieser Methode ist, dass sie schneller ist, vorausgesetzt, die Fol-ge (xn) ist geeignet und die Funktion f ist nicht zu irregulär. Außerdemgibt es ein deterministisches Abbruchkriterium der Berechnung, und zwardie Koksma-Hlawka ungleichung.

3.5 Vergleich der Monte Carlo und Quasi-MonteCarlo Methoden

Die Benutzung der pseudo-zufälligen Zahlen (damit sind die Zahlen gemeint,die mit Zufalsgeneratoren erzeugt werden) für Berechnung von Erwartungs-werten ist eine ziemlich gängige Methode (Monte Carlo). Es funktioniertaber mit einer Durchschnittlichen Fehlerfreiheit. Exaktheit des Ergebnissesist aber ein sehr wichtiger Aspekt, wenn beispielsweise Ingenieure die wahr-scheinlichkeitstheoretischen Modelle benutzen.Wir müssen trotzdem e inige Punkte im Kopf behalten:

1. Für keine Folge, die explizit angegeben werden kann, kann bewiesenwerden, dass sie wirklich zufälligen Charakter haben. Es gibt tieflie-gende theoretische Gründe dafür. (unter anderem von Donald Knuth,in seinem Lebenswerk „The Art of Computer Programming“ beschrie-ben).

2. Wir können nicht verlangen, dass jede Eigenschaft aus der Mengentheo-rie von der Folge erfüllt ist, die wir konstruiert haben, aber wenigstensdie wichtigsten Eigenschaften, wie Eigenschaft einer Normalen Folgeim Sinne von Borel (wenn die Teilfolgen gleicher Länge sich mit glei-cher Frequenz wiederholen), oder χ Quadrat Regel (zeigt, wie gut eine

10

statistisches Modell (unsere Folge) der erwartete Folge entspricht). Al-lerdings gibt es viele Folgen, die diese Eigenschaften erfüllen, aber inwirklichkeit keine guten pseudo-zufälligen Folgen sind.

3. In der Praxis bekommt man die pseudo-zufälligen Folgen auf folgendeWeise:Familien von Zufallsgeneratoren werden durch statistische Tests über-prüft und die erzeugten Folgen, die nicht so gut sind, werden eliminiert.Es gibt eine breite Vielfalt an Tests. (Beschrieben unter anderem vonD.Knuth). Dadurch erhält man Folgen, die als brauchbar angesehenwerden können. Die meisten von ihnen sind endlich (z.B. periodischeFolgen), aber es ist auch möglich unendliche Folgen zu konstruieren.Die gleichverteilten Folgen, die in der quasi-Monte Carlo Methode ver-wendet werden, sind anders in ihrer Natur als die pseudozufälligen Fol-gen in der Monte Carlo Methode. Über diese gleichverteilten Folgenkann man sagen, dass ihre Eigenschaft der gleichverteilung bleibt be-stehen.

4. Was man noch bei beiden Methoden vergleichen kann, sind die Ge-schwindigkeiten der Berechnung. Dieser Vergleich wurde von Sarkarund Prasad (USA) 1987 gemacht. Die beiden Mathematiker haben ge-zeigt, dass für alle zur Zeit bekannten Folgen kleiner Diskrepanz, die aufeinem s-dimensionalen Intervall [0, 1]s definiert sind, wächst die Anzahlder benötigten Simulationen in Abhängigkeit von s sehr schnell, wennsich die Genauigkeit der quasi-Monte Carlo Methode von der Genau-igkeit der Monte Carlo Methode nicht unterscheiden soll. Daher nutztman die quasi-Monte Carlo Methode eher für Probleme kleiner Dimen-sion.

Literaturverzeichnis

[1] Hozumi Morohosi, Masanori F.: A Practical Approach to the ErrorEstimation of Quasi-Monte Carlo Integrations. (1998). – URL http://www.keisu.t.u-tokyo.ac.jp/Research/METR/1998/METR98-10.pdf

[2] Kainhofer, Reinhold: Quasi-Monte Carlo Methoden. (2005).– URL http://reinhold.kainhofer.com/Papers/Kainhofer_Talk_QMC_Wissenswertes_20-06-2005.pdf

[3] N. Bouleau, D. L.: Numerical Methods for Stochastic Processes. NewYork : Wiley, 1994. – Kapitel 2C S

11