Een parallelle multilevel Monte-Carlo-methode voor de simulatie van stochastische partiële differentiaalvergelijkingen

Een parallelle multilevel Monte-Carlo-methode voor de simulatie van stochastische partiële differentiaalvergelijkingen Pieterjan Robbe

Thesis voorgedragen tot het behalen van de graad van Master of Science in de ingenieurswetenschappen: wiskundige ingenieurstechnieken Promotoren: Prof. Dr. ir. S. Vandewalle Prof. Dr. ir. D. Nuyens Assessoren: Prof. Dr. ir. G. Samaey Prof. Dr. R. Vandebril Begeleider: G. Suryanarayana

Academiejaar 2014 – 2015

c Copyright KU Leuven

Zonder voorafgaande schriftelijke toestemming van zowel de promotoren als de auteur is overnemen, kopiëren, gebruiken of realiseren van deze uitgave of gedeelten ervan verboden. Voor aanvragen tot of informatie i.v.m. het overnemen en/of gebruik en/of realisatie van gedeelten uit deze publicatie, wend u tot het Departement Computerwetenschappen, Celestijnenlaan 200A bus 2402, B-3001 Heverlee, +32-16327700 of via e-mail [email protected]. Voorafgaande schriftelijke toestemming van de promotoren is eveneens vereist voor het aanwenden van de in deze masterproef beschreven (originele) methoden, producten, schakelingen en programma’s voor industrieel of commercieel nut en voor de inzending van deze publicatie ter deelname aan wetenschappelijke prijzen of wedstrijden.

Inhoudsopgave Samenvatting

iii

Lijst van Afkortingen en Symbolen

iv

1 Inleiding 1.1 Monte-Carlo-simulatie . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Algoritme voor MLMC-simulatie . . . . . . . . . . . . . . . . . . . . 1.3 Overzicht van de masterproef . . . . . . . . . . . . . . . . . . . . . .

1 2 9 11

2 Multilevel Monte-Carlo voor elliptische stochastische partiële differentiaalvergelijkingen 2.1 Modelprobleem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Samplen uit een stochastisch veld . . . . . . . . . . . . . . . . . . 2.3 Eindige-volumemethode . . . . . . . . . . . . . . . . . . . . . . . 2.4 Resultaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Besluit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . .

. . . . .

13 13 14 18 22 27

3 Analyse van de MLMC-methode voor het modelprobleem 3.1 Parallellisatie . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Parameteranalyse van het modelprobleem . . . . . . . . . . . 3.3 Levelafhankelijke schatters . . . . . . . . . . . . . . . . . . . . 3.4 Besluit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

29 29 32 39 43

4 Multilevel quasi-Monte-Carlo-methoden 4.1 Veralgemeende situering van het modelprobleem 4.2 Oplossingsstrategie . . . . . . . . . . . . . . . . . 4.3 Roosterregels en functieruimten . . . . . . . . . . 4.4 De MLQMC-schatter . . . . . . . . . . . . . . . . 4.5 Het MLQMC-algoritme . . . . . . . . . . . . . . 4.6 Theoretische rekenkost van de MLQMC-schatter 4.7 Resultaten . . . . . . . . . . . . . . . . . . . . . . 4.8 Besluit . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

45 45 46 47 51 54 55 56 60

5 Multilevel Monte-Carlo met meerdere eindfuncties 5.1 Uitbreiding naar meerdere eindfuncties . . . . . . . . . . . . . . . . . 5.2 Numerieke resultaten . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Besluit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61 61 62 66

6 Conclusies en verder werk

69

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

i

Inhoudsopgave A Analytische uitdrukkingen voor eigenwaarden en eigenfuncties in de KL-expansie

73

B De poster

77

C Het artikel

79

Bibliografie

87

ii

Samenvatting In een wiskundig model voor een bepaald probleem of proces zijn de parameters vaak niet exact gekend. In een model voor een wagen bijvoorbeeld, dat wordt gebruikt tijdens een botsproef of crash test, is er een onzekerheid op de geometrie van de wagen door kleine afwijkingen in het productieproces. Enkele voorbeelden zijn het gewicht van de motor, de voorspanning van de bouten, zwakke plekken in de carrosserie. . . Omdat de inputs van het model onzeker zijn zal ook de output onderhevig zijn aan bepaalde onzekerheid. Het domein binnen de ingenieurswetenschap dat zich bezig houdt met het kwantificeren, karakteriseren en reduceren van deze onzekerheid noemt met onzekerheidskwantificatie (uncertainty quantification). In deze masterproef bestuderen we modellen in de vorm van partiële differentiaalvergelijkingen of PDEs. Deze modellen worden wiskundig opgelost door de geometrie te benaderen door een rooster van discrete punten. We onderzoeken en implementeren multilevel Monte-Carlo (MLMC), een recente methode voor het kwantificeren van de onzekerheid in deze PDEs. In de klassieke Monte-Carlo-methode wordt verschillende keren een oplossing berekend voor verschillende realisaties van de onzekere modelparameters. Voor elke realisatie wordt dan een deterministisch probleem (zonder onzekerheid) opgelost. Nadien worden de vele deterministische oplossingen gecombineerd tot een stochastische beschrijving van de oplossing. Het nieuwe MLMC-algoritme combineert de klassieke Monte-Carlo-methode met een sequentie van roosters met verschillende nauwkeurigheid. We tonen aan dat deze nieuwe methode een significante verbetering is ten opzichte van standaard MonteCarlo-technieken. We onderzoeken twee manieren waarop de methode nog verder kan worden versneld. Vooreerst implementeren we het algoritme op een parallelle computer. De vele deterministische problemen kunnen immers allen onafhankelijk van elkaar worden opgelost. Daarnaast onderzoeken we hoe de multilevel Monte-Carlo-methode kan worden gecombineerd met zogenaamde quasi-Monte-Carlo-punten. We tonen aan dat het voor bepaalde problemen mogelijk is om een rekenkost te bereiken die omgekeerd evenredig is met de gevraagde tolerantie op de fout.

iii

Lijst van Afkortingen en Symbolen Symbolen E V a.b a&b ahb Q M` = (m0 2` )d d e( · ) N (N` ) Y` = Q` − Q`−1 C (C` ) α β γ k(x; ω) λ σ2 p(x, ω) C(x, y) mKL PN frac( · ) ∆ q Cd iv

verwachte waarde-operator variantie-operator a/b < constante a/b > constante a/b = constante de quantity of interest, een functionaal Ω → R levelstructuur (aantal roosterpunten op level `) de fysische dimensie van het probleem root-mean-suare-error (RMSE) van een schatter gevraagde tolerantie op de RMSE aantal samples (op level `) verschil van de quantity of interest op level ` en ` − 1 rekenkost (op level `) orde van convergentie van E[|Q` − Q`−1 |] orde van convergentie van V[Q` − Q`−1 ] toename van de rekenkost met M` doorlatendheidstensor (stochastisch veld) correlatielengte variantie druk (stochastisch veld) covariantiefunctie aantal termen in de KL-expansie puntenset met N punten fractioneel deel verschuiving van een puntenset aantal verschuivingen (shifts) multiplicatiefactor voor het aantal samples in MLQMC

Afkortingen MC MLMC QMC MLQMC ODE PDE SPDE SDE RMSE MSE KL-expansie RKHS

Monte-Carlo multilevel Monte-Carlo quasi-Monte-Carlo multilevel quasi-Monte-Carlo gewone differentiaalvergelijking partiële differentiaalvergelijking stochastische partiële differentiaalvergelijking stochastische differentiaalvergelijking root-mean-suare-error mean-square-error, RM SE 2 Karhunen–Loève-expansie Hilbertruimte met reproducerende kern

v

Hoofdstuk 1

Inleiding Bij het opstellen van een wiskundig model voor een proces of probleem in de ingenieurswetenschap en de wetenschap in het algemeen zijn de parameters vaak niet exact gekend [10]. Onzekerheidskwantificatie uncertainty quantification houdt zich bezig met het karakteriseren en reduceren van deze onzekerheid in de modelparameters. Wanneer de input van een model onzeker is, zal immers ook de output onderhevig zijn aan een bepaalde onzekerheid. Bronnen voor deze onzekerheid1 kunnen vaag gespecificeerde of willekeurig variërende parameters zijn. Ook onzekerheid op de geometrie van het probleem, de beginvoorwaarden, materiaaleigenschappen e.d. kunnen onzeker zijn [18]. In deze masterproef bespreken we hoe de multilevel Monte-Carlo-methode (MLMC) kan worden gebruikt voor het simuleren van modellen beschreven door partiële differentiaalvergelijkingen. Monte-Carlo-methoden zijn ontstaan in de loop van de nucleaire proeven tijdens het Manhattanproject in 1943 [33]. De ontwikkelaars van de methode, Stan Ulam en John von Neumann, suggereren dat de gebruikte computerexperimenten kunnen worden versneld door de botsingen tussen de gesimuleerde deeltjes te benaderen met een statistisch equivalent. Net als het hele Manhattanproject bleef ook de nieuwe methode geheim. Door Nicholas Metropolis werd het pseudoniem Monte Carlo gebruikt, naar het gelijknamige casino in Monaco, waar de oom van Stan Ulam een fervent gokker was. De kern van een Monte-Carlo-simulatie is het voorstellen van de te berekenen grootheid als de verwachte waarde van een stochastisch proces, en vervolgens deze grootheid te benaderen aan de hand van een groot aantal realisaties van dat stochastische proces. In de jaren daarna won het Monte-Carlo-algoritme aan populariteit, voornamelijk in de financiële wereld voor het simuleren van aandelenkoersen en optieprijzen. Het grootste nadeel aan de methode is de trage convergentie: er zijn zeer veel realisaties van het stochastisch proces vereist om de oplossing voldoende nauwkeurig te berekenen. Verschillende verbeteringen werden aangebracht om de convergentie van het algoritme te versnellen door middel van variantiereductietechnieken. Daarnaast werden heel wat afgeleide technieken ontwikkeld zoals Markov 1

Fysisch maakt men vaak een onderscheid tussen inherente onzekerheid (toeval) of epistemische onzekerheid (gebrek aan kennis) [6].

1

1. Inleiding chain Monte Carlo (MCMC), Metropolis-Hastings. . . In 2008 werd door M. Giles het multilevel Monte-Carlo-algoritme geformuleerd als een mogelijke variantiereductietechniek [17]. De methode gebruikt een hiërarchie van geneste roosters. De verwachte waarde van het stochastisch proces op het fijnste rooster kan dan worden geschreven als de verwachte waarde van dat stochastisch proces op een ruw rooster, aangevuld met een reeks correctietermen. Omdat de variantie van de correctietermen in de som klein zijn heeft men minder én goedkopere realisaties nodig dan de standaard Monte-Carlo-schatter. Dit multilevelidee werd eerder al voorgesteld door A. Kebaier [24] en onafhankelijk ook door S. Heinrich [22] in het kader van parametrische integratie. In een eerste deel van deze masterproef beschouwen we de theoretische fundamenten van de MLMC-methode. We bespreken standaard Monte-Carlo-technieken en multilevel Monte-Carlo-technieken in hun meest algemene vorm, stellen een algoritme op en proberen een schatting te maken van de rekenkost.

1.1

Monte-Carlo-simulatie

Bij Monte-Carlo-simulatie (MC) is men geïnteresseerd in de verwachte waarde van een functionaal Q = G(X) van X, een multivariate toevalsvariabele of toevalsvector in de oneindigdimensionale kansruimte (Ω, F, P). Ω is de uitkomstenruimte, de verzameling van alle mogelijke uitkomsten van een kansexperiment; F is een gebeurtenis als deelruimte van Ω en P is een kansmaat met P(Ω) = 1. Deze functionaal Q : Ω → R, ook wel de eindfunctie of quantity of interest genoemd, kan een lineaire of nietlineaire functie zijn van X (zie §3.2.3). Wanneer X berekend wordt als oplossing van een vergelijking is X vaak niet exact beschikbaar. De oplossing X kan dan worden benaderd door een eventueel discrete toevalsvariabele XM . In het geval van een PDE bijvoorbeeld, komt M overeen met de grootte van de ruimtelijke discretisatie. De functionaal wordt dan QM = GM (XM ), waarbij we aannemen dat de verwachte waarde van het verschil tussen Q en QM naar nul convergeert voor M → ∞. Doorgaans is men dan geïnteresseerd in het vinden van de parameter α zodat |E[QM − Q]| . M −α , (1.1) en in het opstellen van een methode die deze parameter zo groot mogelijk maakt. Met de notatie “a . b” wordt bedoeld dat ab < c met c een zekere constante die niet afhangt van andere parameters en waarvan abstractie wordt gemaakt. Een analoge definitie geldt voor a & b en a h b. De orde van convergentie in (1.1) is dan α. Voor het voorbeeld van de PDE wensen we dat de verwachte waarde van het verschil tussen een functionaal van de oplossing op het rooster en de functionaal van de oneindigdimensionale oplossing naar nul gaat naarmate de ruimtelijke discretisatiestap verkleint, en dat deze waarde naar boven begrensd is door c · M −α . We zijn geïnteresseerd in het schatten van E[Q]. Gegeven M dan hzijn we op zoek i ˆ M voor E[QM ] waarvoor geldt dat E Q ˆ M − E[Q] = naar benaderingen of schatters Q 0. We kwantificeren de nauwkeurigheid van deze benaderingen aan de hand van de 2

1.1. Monte-Carlo-simulatie root-mean-square error (RMSE) ˆ M − E[Q])2 ˆ M , E (Q e Q

h

i1/2

.

In wat volgt bespreken we Monte-Carlo-simulatie en multilevel Monte-Carlo-simulatie, ˆ M voor E[QM ]. die enkel verschillen in de gebruikte schatter Q

1.1.1

Het standaard Monte-Carlo-algoritme

De standaard Monte-Carlo-schatter voor E[QM ] is N X ˆ MC , 1 QM (ωn ), Q M,N N n=1

(1.2)

zie [5]. Het is de gemiddelde waarde van een steekproef van N onafhankelijke trekkingen (samples) QM (ωn ), met ωn een realisatie van ω, een stochastische parameter, die aangeeft dat de functionaal QM een stochastische variabele is. Er zijn twee oorzaken van fouten in de schatter (1.2): de benadering van Q door de eindigdimensionale functionaal QM , gerelateerd aan de ruimtelijke discretisatie, en de stochastische fout door het benaderen van de verwachte waarde door het steekproefgemiddelde. De bijdrage van deze fouten in de RMSE kan men vinden door de mean-square-error (MSE = RMSE2 ) uit te schrijven: ˆ MC e Q M,N

2

2

=E

ˆ MC − E[Q] Q M,N

=E

ˆ MC − E[Q] ˆ MC − E Q ˆ MC + E Q Q M,N M,N M,N

=E

ˆ MC ˆ MC Q M,N − E QM,N

h

i

h

h

=E

ˆ MC − E Q ˆ MC Q M,N M,N h

i

2 2

ˆ MC + E Q M,N − E[Q] h

ˆ MC ˆ MC − E Q 2 Q M,N M,N

i2

h

i2

i

+

ˆ MC − E[Q] E Q M,N

i h

i

i

2

ˆ MC − E[Q] + E Q M,N h

i

ˆ MC + (E Q ˆ MC − E[Q])2 . =V Q M,N M,N h

i

h

i

ˆ MC een zuivere schatter is, In de derde gelijkheid valt het kruisproduct weg omdat Q M,N h h ii MC MC ˆ ˆ i.e. E Q −E Q = 0. Omdat verder geldt dat M,N

E

h

M,N

ˆ MC Q M,N

en V is

h

ˆ MC Q M,N

i

i

N N 1 X 1 X =E QM (ωn ) = E [QM (ωn )] = E [QM ] N n=1 N n=1

"

#

N N 1 X 1 X 1 =V QM (ωn ) = 2 V [QM (ωn )] = V [QM ] N n=1 N n=1 N

"

ˆ MC e Q M,N

#

2

= N −1 V [QM ] + (E [QM − Q])2 .

(1.3) 3

1. Inleiding De eerste term is de variantie van de Monte-Carlo-schatter, die de stochastische fout voorstelt door het benaderen van de exacte waarde van de functionaal door een eindige√som, en afneemt met N −1 . De RMSE van een Monte-Carlo-schatter gaat dus met 1/ N , een bekend resultaat [5]. De tweede term is afkomstig van de ruimtelijke discretisatiefout, d.w.z. het verschil tussen QM en Q. ˆ MC kleiner Een voldoende voorwaarde voor het bereiken van een RMSE e Q M,N dan is dat beide termen in het rechterlid van (1.3) kleiner zijn dan 2 /2. Voor de stochastische fout dient men een voldoende groot aantal samples N te nemen. Voor de tweede term moet de ruimtelijke discretisatie fijn genoeg zijn, zodat E [QM − Q] = ˆ MC O(). In de veronderstelling dat V [QM ] onafhankelijk is van M , is de schatter Q M,N een voldoende nauwkeurige benadering voor E[Q] als N −1 V [QM ] & 2 /2, of nog, N & −2 en, gebruik makende van (1.1), als M & −1/α . Veronderstel dat de kost voor het berekenen van één sample QM (ω) naar boven begrensd is als volgt C (QM (ω)) . M γ ,

voor γ > 0. De kost van de Monte-Carlo-schatter met N samples is dan ˆ MC . N M γ . C Q M,N

De kost voor het bereiken van een RMSE van , ook wel de -kost genoemd, voldoet dan aan de volgende vergelijking: ˆ MC . −2−γ/α , C Q M,N

(1.4)

zie [8], immers N & −2 en M & −1/α .

1.1.2

De multilevel Monte-Carlo-methode

Het basisidee van multilevel Monte-Carlo-simulatie (MLMC) [17] is eenvoudig: neem samples niet enkel uit één benadering QM voor Q maar uit meerdere benaderingen QM` , ` = 0 . . . L. Beschouw de geometrische reeks {M` ∈ N : ` = 0 . . . L}

(1.5)

waarbij M` = s · M`−1 voor ` ≥ 1 en met s een natuurlijk getal groter dan één: M0 < M1 < . . . < ML . De M` -waarden zullen doorgaans overeenstemmen met het aantal roosterpunten in het discretisatierooster van een PDE. Deze ` noemt men dan levels. Ze corresponderen met discretisatieroosters met een verschillende fijnheid van discretisatiegrootte. Het multilevelidee bestaat erin om de verwachte waarde E [QM` ] niet rechtstreeks te schatten, maar eerder als de verwachte waarde op het vorige level ` − 1 plus de verwachte waarde van een zekere correctieterm Y` , QM` − QM`−1 , of i h nog, E [QM` ] = E QM`−1 + E [Y` ]. Dit idee kan recursief worden doorgevoerd, als volgt. Startend van E [QM ] = E [QM0 ] − 4

L X `=1

E [QM` ] + E[QM`−1 ] ,

1.1. Monte-Carlo-simulatie en wegens lineariteit van de verwachte waarde, geldt dat E [QM ] = E [QM0 ] +

L X `=1

E QM` − QM`−1 = i

h

L X

E [Y` ] .

(1.6)

`=1

met Y0 = QM0 . Met andere woorden, de verwachte waarde van de gezochte functionaal op het fijnste rooster (` = L) is gelijk aan de verwachte waarde op het ruwste rooster (` = 0) plus een som van correctietermen die het verschil zijn in verwachte waarde tussen de opeenvolgende levels. Men schat dus niet de verwachte waarde E [QM` ] op elk level, maar wel de verwachte waarde van opeenvolgende verschillen Y` (ω). Merk op dat we in de vorige notatie nogmaals de expliciete afhankelijkheid van de stochastische parameter ω aangeven. We kiezen voor deze schatter eveneens de Monte-Carlo-schatter (1.2) N` 1 X MC (ω ) (ω ) − Q Q Yˆ`,N , n n M M `−1 ` ` N` n=1

(1.7)

met N` het aantal samples op elk level. De multilevel Monte-Carlo-schatter (verder MLMC-schatter genoemd) kan dan gedefinieerd worden als volgt: ˆ MLMC , Q M,{N` }

L X

MC Yˆ`,N . `

(1.8)

`=0

Voor Y` kan eender welke zuivere schatter (unbiased estimator) worden gebruikt, zoals een quasi-Monte-Carlo-schatter [5, 21] (zie Hoofdstuk 4) of een meer klassieke kwadratuur- of cubatuurregel. Belangrijk is wel dat de schatting (1.7) gebeurt op basis van dezelfde stochastische parameter ω ∈ Ω op beide levels M` en M`−1 . Eenzelfde foutenanalyse als voor de Monte-Carlo-schatter levert ˆ MLMC e Q M,{N` }

2

=E =

L X `=0

ˆ MLMC − E[Q] Q M,{N` }

2

N`−1 V[Y` ] + (E[QM − Q])2 ,

(1.9)

immers, waarden E[Y` ] worden onafhankelijk van elkaar geschat en dus h de verwachte i PL MLMC ˆ is V QM,{N` } = `=0 N`−1 V[Y` ]. De fout op de MSE (1.9) bestaat dus opnieuw uit twee delen: een stochastische fout voorgesteld door de variantie van de MLMCschatter en een discretisatiefout. Om een RMSE< te bekomen moeten beide termen opnieuw kleiner zijn dan 2 /2. Merk op dat de discretisatiefout dezelfde is als voor de Monte-Carlo-schatter. Men kiest dus opnieuw M = ML & −1/α . Een voorwaarde voor N` volgt in §1.1.3. DeMLMC-schatter is pas voordelig wanneer het goedkoper is om te voldoen MLMC MC ˆ ˆ aan e QM,{N` } < , zie (1.9), dan aan e QM,N < , zie (1.3). Dit is het geval wanneer men minder en/of goedkopere samples nodig heeft om de variantie van de schatter onder de 2 /2-grens te houden. Dit blijkt inderdaad zo te zijn voor de MLMC-schatter (1.8) als twee voorwaarden zijn voldaan [8]: 5

1. Inleiding 1. De variantie van Y` gaat naar 0 als ` → ∞. Er zijn dus steeds minder samples op elk level ` nodig om E [Y` ] te schatten naarmate ` groter wordt. Dit is de reden voor de voorwaarde dat de schattingen voor Q` (ω) en Q`−1 (ω) moeten gebaseerd zijn op dezelfde stochastische parameter ω ∈ Ω. 2. Voor voldoende kleine Mi , met i < `, zijn samples van QMi veel goedkoper dan samples van QM` .

1.1.3

Rekenkost van de MLMC-schatter

De totale -rekenkost voor de MLMC-schatter wordt gegeven door ˆ MLMC = C Q M,{N` }

L X `=0

(1.10)

N` C` ,

waarbij C` de kost is voor het nemen van één sample van Y` (ω). Om te bepalen hoeveel samples N` er op elk level ` nodig zijn, kan men volgend optimalisatieprobleem beschouwen [18]. Men wenst de totale rekenkost van de MLMC-schatter te minimaliseren waarbij de variantie van de schatter binnen de gevraagde tolerantie blijft. Het optimalisatieprobleem kan als volgt worden geformuleerd: min N`

s.t.

X `

X

N` C` N`−1 V[Y` ]

`

2 = 2

(1.11) .

Merk op dat het aantal samples op elk level als een continue variabele wordt gezien. Bovenstaande formulering kan worden omgezet in een probleem zonder beperkingen via de methode van de Lagrangevermenigvuldigers. Hiertoe definieert men de zogenaamde Lagrangiaan, L(N` , λ1 ) =

X `

N` C` + λ1 ·

X

N`−1 V[Y` ]

`

2 − 2

!

(1.12)

waarbij λ1 een nieuwe variabele in het probleem is. De oplossing van het oorspronkelijk optimalisatieprobleem (1.11) is een kritisch punt van de Lagrangiaan (1.12). Door het toepassen van de eerste-orde optimaliteitsvoorwaarde (de partiële afgeleide van de Lagrangiaan (1.12) naar N` moet nul zijn) verkrijgt men dat ∂L V[Y` ] = C` − λ 1 =0 ∂N` N`2 of

6

⇔ C` N`2 = λ1 V[Y` ]

N` =

p

h

q

q

λ1 V [Y` ] C`−1 V [Y` ] C`−1

(1.13)

1.1. Monte-Carlo-simulatie voor elk level `. Uit de tweede partiële afgeleide van de Lagrangiaan vinden we de gelijkheidsbeperking terug: X 2 ∂L −1 = N` V[Y` ] − = 0. ∂λ1 2 `

Combinatie met (1.13) levert volgende vergelijking: X

N`−1 V[Y` ] =

`

s X

⇔

`

C` 2 V[Y` ] = λ1 V[Y` ] 2 s

X

⇔

2 2

`

C` V[Y` ] 2 = λ1 2 λ1 =

p

⇔

2 Xq C` V[Y` ]. 2 `

(1.14)

Merk op dat (1.13) hetzelfde verband geeft als wanneer men de variantie van de schatter voor E [Y` ] minimaliseert, waarbij de rekenkost op elk level constant blijft. De Lagrangiaan wordt in dit geval L(N` , λ2 ) = V Yˆ` + λ2 · (N` C` − cte) h

=

i

X 1 `

N`

V [Y` ] + λ2 · (N` C` − cte) ,

en het kritisch punt vinden we als ∂L 1 = − 2 V[Y` ] + λ2 C` = 0 ∂N` N` ⇔ λ2 C` N`2 = V[Y` ]

of

1 q N` = √ V [Y` ] C`−1 . λ2

Met de relaties (1.13) en (1.14) wordt de totale rekenkost voor de MLMCschatter (1.10) C

ˆ MLMC = Q M,{N` }

L X `=0

=

p

N` C`

λ1

L q X

V[Y` ]C`

`=0

2 = 2

L q X

V[Y` ]C`

!2

.

`=0

7

1. Inleiding Wanneer de variantie V[Y` ] sneller daalt met ` dan de rekenkost C` stijgt dan ligt de dominante rekenkost op het ruwste rooster ` = 0. Dat betekent dat de rekenkost van het MLMC-algoritme proportioneel zal zijn met V[Y0 ]C0 . De rekenkost voor gewone Monte-Carlo is evenredig met V[Y0 ]CL . Er is immers slechts één (fijnste) rooster waarop de samples worden genomen. De winst aan rekenkost van multilevel MonteCarlo ten opzichte van de standaard Monte-Carlo-schatter is dan proportioneel met (C0 /CL ) h (M0 /ML )γ h γ/α . Omgekeerd, wanneer de variantie V[Y` ] sneller stijgt met ` dan de rekenkost C` daalt, dan ligt de dominante rekenkost op het fijnste level ` = L. De rekenkost van MLMC-schatter is proportioneel met V[YL ]CL . De winst in rekenkost zal dan bepaald worden door de verhouding van de varianties V[YL ]/V[Y0 ]. De resultaten voor de -rekenkost van de MLMC-schatter worden samengevat in onderstaand Theorema [1]. Voor het bewijs verwijzen we naar [8, 17]. Theorema 1. Veronderstel dat er positieve constanten α, β en γ bestaan met α ≥ min(β, γ) en dat volgende drie voorwaarden zijn voldaan: (i) |E [QM` − Q]| . M`−α ,

(ii) V[Y` ] . M`−β en (iii) C` . M`γ .

Dan geldt voor elke < exp(−1) dat er een getal L en een rij {N` }L `=0 bestaan zodat de MSE 2 2 MLMC MLMC ˆ ˆ e QM,{N` } , E QM,{N` } − E[Q] < 2 en dat ˆ MLMC . C Q M,{N` }

 −2   

−2 (log )2

  −2−(γ−β)/α

als β > γ, als β = γ, als β < γ.

(1.15)

Wanneer de probleemafhankelijke parameters α, β en γ gekend zijn, dan kan men voor een gegeven nauwkeurigheid een a priori schatting maken van de kost van het algoritme. Het eerste geval, β > γ in Theorema 1, komt overeen met een variantie die sneller afneemt dan de rekenkost stijgt. De dominante rekenkost ligt op de ruwe roosters waar C` = O(1). Er zijn O(−2 ) samples nodig om de gevraagde nauwkeurigheid te bereiken. Voor het tweede geval, β < γ in Theorema 1, neemt de rekenkost sneller toe met kleiner wordende tolerantie dan de variantie daalt. De dominante rekenkost ligt dan op de fijnste roosters. Omwille van de eerste voorwaarde in Theorema 1 is 2−αL = O() en dus is CL = O(−γ/α ). Als β = 2γ dan is de totale rekenkost O(CL ) en worden O(1) samples genomen op level L, het fijnste rooster. Voor het limietgeval β → 0 benadert de kost van de MLMC-schatter de kost van de MC-schatter. Voor het grensgeval β = γ is de rekenkost en de bijdragen aan de totale variantie gelijk verdeeld over alle levels. In [17] wordt (1.6) een “telescopische som” genoemd. Men zoomt als het ware telkens in op verschillen tussen levels die dan uiteindelijk ineenschuiven om het echte 8

1.2. Algoritme voor MLMC-simulatie resultaat te vormen. Het idee om een geometrische sequentie van stapgroottes te gebruiken is afkomstig uit de praktijk van de multigrid methoden [3] voor het oplossen van partiële differentiaalvergelijkingen. Op het fijnste rooster is de discretisatiefout het kleinst, maar de kost voor een iteratie met bijvoorbeeld Gauss-Seidel is groot. Op het ruwste rooster is de discretisatiefout veel groter, maar de rekenkost ook veel kleiner. Het multilevelalgoritme (1.6) kan ook geïnterpreteerd worden als recursieve controlevariabelenstrategie: elke stochastische grootheid op level ` wordt gecontroleerd door zichzelf op een lager level ` − 1. Merk op dat een geometrische sequentie van levels M` = {1, 2, 4, 8 . . .} niet steeds optimaal is [18]. Hoewel de resultaten die in Hoofdstuk 2 zullen besproken worden aangeven dat MLMC-simulatie erg effectief is, zou het kunnen dat een andere sequentie toch nog betere resultaten oplevert.

1.2

Algoritme voor MLMC-simulatie

Om tot een werkend algoritme te komen moeten nog enkele componenten verder worden aangevuld. Zo moeten heel wat parameters nog een waarde krijgen: M0 (het aantal roosterpunten op het ruwste rooster), s (de geometrische factor in de levelstructuur) en L (het aantal levels). In wat volgt bespreken we ook een praktische formule voor het aantal samples N` op elk level en voor het stopcriterium. De overige parameters moeten oordeelkundig worden gekozen. • Vaak neemt men de geometrische factor s = 2 omdat het een natuurlijke verfijning is van een discretisatierooster van een PDE. Voor bepaalde SDEs toont Giles [17] aan dat de optimale waarde voor s gelijk is aan 7. Dat betekent dat de roosterafstand in level ` 7 keer kleiner is dan op level ` − 1. Giles raadt dan aan om s ≤ 4 te gebruiken als vuistregel. Op die manier heeft men nog steeds de voordelen van een grote waarde voor s en zijn een voldoende aantal levels mogelijk om de systematische fout (bias) te schatten. • M0 hangt af van de maximaal toelaatbare stapgrootte in de ruimtelijke discretisatie van de PDE in de gegeven toepassing. • Het aantal levels L wordt bepaald door het fijnste rooster dat computationeel aanvaardbaar is in de PDE toepassing. Uiteraard hangt deze waarde af van M0 en s. Hoe meer levels men kan toevoegen, hoe groter de winst aan rekenkost door het gebruik van een multilevel Monte-Carlo-schatter. Tot slot worden alle elementen samengebracht in een algoritme voor MLMC-simulatie.

1.2.1

Een praktische formule voor N`

Om de RMSE < te houden moet de eerste term (1.9) kleiner zijn dan 2 /2. Het √ inp optimale aantal samples op elk level ` is N` = λ1 V [Y` ] /C` , zie (1.13). Invullen 9

1. Inleiding van (1.14) voor

√

λ1 geeft L q X 2q N` = 2 V[Y` ]/C` V[Yi ]Ci , i=0

&

'

(1.16)

zie [17]. Het optimale aantal samples voor elk level ` = 0 . . . L kan dus pas worden berekend als de variantie V[Y` ] gekend is. In de praktijk neemt men N 0 initiële samples van Y` en schat men hieruit de variantie op level `. Daarna kan formule (1.16) worden geëvalueerd en kan men het gewenste aantal samples N` op elk level ` bepalen.

1.2.2

Stopcriterium

De optimale waarde van het aantal levels L kan “on the fly” worden berekend uit de steekproefgemiddeldes en -varianties. Veronderstel dat er een getal M 0 bestaat zodanig dat de convergentie van het linkerlid in (1.1) monotoon is voor elke M ≥ M 0 . Als |E[Y` ]| = |E[Q` − Q`−1 ]| ∝ 2−α` (zie Theorema 1) dan is de overblijvende fout E[Q − QL ] = ≈ ≈

∞ X

`=L+1 ∞ X

E[Q` − Q`−1 ] 2−α`

`=L+1 ∞ X −αL −αk

2

2

k=1

≈ 2−αL

∞ X

2−αk

k=1

≈ E[QL − QL−1 ]

∞ X

2−αk

k=1

1 ≈ E[QL − QL−1 ] −1 1 − 2−α

1 ≈ E[QL − QL−1 ] α 2 −1

!

!

waarbij de somformule voor een oneindige reeks werd gebruikt. Er is aan de gevraagde √ tolerantie op de RMSE voldaan als |E[Q − QL ]| < / 2, en dus volgt dat |E[YL ]| = |E[QL − QL−1 ]| < √ (2α − 1) . 2

(1.17)

Door extrapolatie kan men ook het vorige punt |E[YL−1 ]| in rekening brengen. Zo verkrijgt men een meer robuust stopcriterium, dat later in de numerieke experimenten wordt gebruikt: max{2−α |YˆL−1 |, |YˆL |} < √ (2α − 1) , (1.18) 2 zie [17, 18]. 10

1.3. Overzicht van de masterproef

1.2.3

Algoritme: MLMC-simulatie

Alle voorgaande bemerkingen samengenomen kan een algoritme voor MLMC-simulatie worden opgesteld [8]: 1 2 3 4

5 6 7 8 9 10

Start met L = −1, geconvergeerd = false; repeat L ← L + 1; Schat V[YL ] door de steekproefvariantie voor een initieel aantal samples N 0; Bereken het optimale aantal samples op elk level ` = 0 . . . L met (1.16); Neem extra samples op elk level ` ≤ L volgens de nieuwe N` ; if L ≥ 1 then Test op convergentie met (1.18); end until geconvergeerd;

Stap (6) probeert de variantie van de schatter kleiner dan 2 /2 √ te houden, terwijl stap (8) er voor zorgt dat de resterende bias kleiner is dan / 2. Samen moeten deze twee beperkingen de fout op de RMSE binnen de gevraagde tolerantiegrens houden. Het initieel aantal samples voor de schatting van de variantie kiezen we als een constante, N 0 . Merk op dat het algoritme heuristisch is. Het is immers niet wiskundig gegarandeerd dat het uiteindelijk berekende resultaat een RMSE heeft kleiner dan de gevraagde tolerantie . De meest delicate component van het algoritme is de schatting van de fout E [Q − QL ] via (1.18). Een oordeelkundige keuze voor de discretisatiefijnheid op het fijnste rooster, level L, vergt een grondige kennis van het probleem. Ook bij het oplossen van deterministische PDEs is het schatten van de fout en de constructie van een optimaal (eventueel adaptief) verfijnd rooster een moeilijk probleem. De oplossing van dat probleem is nog volop onderwerp van onderzoek in heel wat onderzoeksgroepen. Dit ligt dan ook buiten het bestek van deze masterproef. In elk geval is een overschatting van het fijnste level L minder erg voor multilevel Monte-Carlo dan voor gewone Monte-Carlo. De MLMC-schatter werkt immers op L verschillende levels, en een te grote L resulteert enkel in een hogere kost voor het fijnste rooster. De gewone MC-schatter zou enkel samples nemen op dat fijne (en dus dure) discretisatierooster.

1.3

Overzicht van de masterproef

Het vervolg van deze tekst is als volgt ingedeeld. In Hoofdstuk 2 bespreken we hoe het multilevel Monte-Carlo-algoritme kan worden toegepast op het modelprobleem van een elliptische stochastische partiële differentiaalvergelijking. Deze vergelijkingen worden bijvoorbeeld gebruikt bij het modelleren van grondwaterstroming doorheen poreuze media. In Hoofdstuk 3 tonen we aan dat het multilevelalgoritme zich uitstekend leent tot parallelle uitvoering, net als standaard Monte-Carlo. Met behulp van de parallelle implementatie van het algoritme kan men de uitvoeringstijd tot 18 11

1. Inleiding keer versnellen. Voor een voldoende uitdagend probleem, dus kleine tolerantie, is de enige beperking het aantal processoren dat men ter beschikking heeft. Verder voeren we een parameteranalyse uit van het modelprobleem en bestuderen we levelafhankelijke schatters. In de recente literatuur probeert men het multilevel Monte-Carlo-algoritme nog verder te versnellen door het te combineren met bestaande variantiereductietechnieken. Enkele voorbeelden zijn MLMC met importance sampling [4] en MLMC met antithetische variabelen [19]. Hoofdstuk 4 bevat de belangrijkste eigen bijdrage. We tonen aan dat men door het combineren van multilevel Monte-Carlo en quasi-MonteCarlo, als variantiereductietechniek, een substantiële vermindering in rekenkost kan bekomen. In een laatste hoofdstuk bekijken we tenslotte hoe het multilevelalgoritme kan worden aangepast voor het berekenen van meerdere eindfuncties in éénzelfde simulatie.

12

Hoofdstuk 2

Multilevel Monte-Carlo voor elliptische stochastische partiële differentiaalvergelijkingen In dit hoofdstuk wordt een veel gebruikt modelprobleem voor de simulatie van stochastische partiële differentiaalvergelijkingen (SPDEs) besproken. Dergelijke vergelijkingen worden bijvoorbeeld gebruikt bij onzekerheidskwantificatie voor grondwaterstroming. De onzekerheid in het probleem is de diffusiecoëfficiënt die kan worden gemodelleerd als een stochastisch veld. Één realisatie van het stochastisch proces behelst dan (i) een sample nemen van het stochastisch veld en (ii) de deterministische partiële differentiaalvergelijking (PDE) oplossen. Voor het eerste deelprobleem bekijken we de Karhunen–Loève-expansie en leiden analytische uitdrukkingen af voor de eigenwaarden en eigenfuncties die in deze expansie nodig zijn. De PDE wordt dan in de daarop volgende stap opgelost met behulp van een eindige-volumemethode. We bestuderen de MLMC-schatter als variantiereductietechniek voor dit modelprobleem en tonen aan de hand van numerieke experimenten in één en twee dimensies de superioriteit ten opzichte van standaard Monte-Carlo-methoden.

2.1

Modelprobleem

Het beschouwde probleem is de stroming van grondwater doorheen een poreus medium [40]. Dergelijke stroming wordt beschreven door de wet van Darcy, genoemd naar de Franse wetenschapper Henry Darcy die in 1856 zijn eerste bevindingen publiceerde. Zijn experimenten toonden aan dat het debiet van een laminaire stroming doorheen een volume recht evenredig is met het verschil in waterhoogte aan de beide einden van dit volume en omgekeerd evenredig met de lengte ervan [15]. Vanwege de lage stroomsnelheden in grondwater is het laminaire karakter van de stroming een geldige aanname. Het model voor grondwaterstroming kan eveneens worden afgeleid als een vereenvoudiging van de algemene Navier-Stokes vergelijkingen. De onzekerheid in het probleem is de doorlatendheid (hydraulic conductivity) k, een maat voor de weerstand die een vloeistof ondervindt bij stroming door een 13

2. Multilevel Monte-Carlo voor elliptische SPDEs poreus medium. Als mogelijke toepassing vermelden we de risicoanalyse bij de opslag van radioactief materiaal (of een andere, gevaarlijke of polluerende stof). Grondwater dat in contact komt met gelekt radioactief afval, kan in contact komen met de mens na transport door het aardoppervlak. De kenmerken van een dergelijk grondwatertransport kunnen moeilijk experimenteel worden opgemeten. Omwille van de grote tijdsschalen in het proces moet men zich behelpen met modellen en simulaties. De klassieke vergelijkingen voor een stationaire éénfase-stroming bestaan uit de wet van Darcy gekoppeld met het divergentievrij snelheidsveld van een onsamendrukbare stroming: q + k∇p = g (2.1) ∇·q =0 in een gebied D ∈ Rd . In de praktijk is de dimensie d natuurlijk gelijk aan 3. Omwille van de hoge rekenkost gebruikt men echter ook vaak tweedimensionale (d = 2) of zelfs ééndimensionale modellen (d = 1). In (2.1) is q de Darcyflux [m/s], ∇p de drukgradiënt [Pa/m] en k = κ/µ de doorlatendheidstensor. κ is de permeabiliteit gemeten in [m2 ] of eenheden “darcy” en µ is de viscositeit in [kg/m/s]. Het rechterlid g is een bronterm. Bemerk de analogie met de wet van Ohm in de elektriciteitsleer: de doorlatendheid is omgekeerd evenredig met de weerstand, en de stroming q is het potentiaalverschil, ∇p, gedeeld door die weerstand. Met de aanname dat k een stochastisch veld is (zie hieronder) kan het modelprobleem (2.1) worden herschreven tot een tweede-orde elliptische partiële differentiaalvergelijking − ∇ · (k(x; ω)∇p(x; ω)) = f (x; ω) (2.2) met f (x; ω) = −∇ · g de divergentie van de gekende bronterm. De gezochte oplossing p(x; ω) is eveneens een stochastisch veld op D × Ω. In wat volgt beperken we ons tot D = [0, 1]d . Het nemen van een sample uit een SPDE bestaat altijd uit twee delen:

1. Neem een sample van de stochastische parameter ω en herleid de SPDE naar een deterministische PDE, en 2. los het deterministisch probleem op. In wat volgt bespreken we de modellering van k als een stochastisch veld (random field) en tonen we aan hoe k kan geschreven worden als een lineaire combinatie van orthogonale basisfuncties. Daarna bespreken we de eindige-volumemethode in één en twee dimensies als een oplossingsstrategie voor een deterministische PDE.

2.2

Samplen uit een stochastisch veld

De enige stochastische parameter in het modelprobleem (2.1) is de doorlatendheid k. Vaak wordt deze gemodelleerd als een stochastisch veld k(x; ω) ∈ D × Ω. Dit is de natuurlijke uitbreiding van een stochastisch proces met de tijd als enige dimensie, 14

2.2. Samplen uit een stochastisch veld 1 0.8

x2

0.6 0.4 0.2 0

0

0.2

0.4

0.6

0.8

1

x1 Figuur 2.1: Een typische realisatie van het lognormaal verdeelde veld k(x; ω) voor σ 2 = 1 en λ = 0.01 en D = [0, 1]2 .

zoals bijvoorbeeld een Brownse beweging, naar het meerdimensionale geval met ruimtelijke dimensies. Het veld wordt gekarakteriseerd door een gemiddelde waarde en een covariantiefunctie C(x, y). Beiden moeten worden geschat uit de data die voorhanden is. Vaak neemt men een lognormale verdeling aan voor k: een veld met scalaire waarden waarvan de logaritme Gaussiaans is. Dit is een niet-triviale distributie die computationeel voordelen biedt maar toch voldoende dicht bij de gemeten data aansluit. De doorlatendheid is immers steeds positief en varieert snel, zelfs over kleine lengte-schalen. Figuur 2.1 toont een typische realisatie van een lognormaal verdeeld veld, in twee dimensies. In de literatuur stelt men volgende vorm voor de covariantiefunctie voorop [8]: kx − ykp C(x, y) , σ exp − λ 2

voor x, y ∈ D.

(2.3)

Hierin is k · kp de klassieke p-norm in Rd , σ 2 de variantie van het stochastisch veld en λ de correlatielengte. Deze laatste bepaalt hoe sterk de waarden van het veld in twee nabijgelegen punten van het veld met elkaar gecorreleerd zijn. Het geeft aan hoe groot de invloedssfeer van een welbepaalde positie is. Hoe groter de variantie en kleiner de correlatielengte, hoe grilliger het veld, en dus hoe moeilijker de op te lossen SPDE (zie Hoofdstuk 3). Voor grondwaterstroming is men typisch geïnteresseerd in toepassingen met λ ≤ diam(D) en σ 2 ≥ 1. De covariantiefunctie is stationair (een verschoven versie levert hetzelfde resultaat op) en isotroop (ze hangt enkel of van de afstand |x − y| tussen de twee punten). Zo’n functie C wordt een radiale basisfunctie (RBF) genoemd. De covariantiefunctie in één dimensie is weergegeven in Figuur 2.2. Er bestaan verschillende methoden voor het nemen van samples uit een stochastisch veld met gegeven covariantiefunctie. De drie belangrijkste zijn wellicht de polynomial chaos expansion [11], de Karhunen–Loève (KL)-expansie [16] en circu15

2. Multilevel Monte-Carlo voor elliptische SPDEs

1

C(x, y)

C(x, y)

1

0.5

1

1

1 0.5

y

0.5

1 0.5

0.5 0 0

y

x

0.5 0 0

x

Figuur 2.2: (a) Exacte covariantie-oppervlak C(x, y), 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 voor een exponentiële covariantie (2.3) in de 1-norm met correlatielengte λ = 1 en variantie σ 2 = 1. (b) 4-terms benadering door middel van een KL-expansie. lant embedding [21]. We gebruiken hier de KL-expansie1 voor het stochastisch veld Z(x; ω) , log k(x; ω). Z(x; ω) is dus een Gaussiaans verdeeld veld. Men wil nu Z schrijven als een oneindigdimensionale lineaire combinatie van orthogonale functies van de vorm Z(x; ω) = Z¯ +

∞ p X

θn fn (x)ξn (ω),

(2.4)

n=1

analoog aan de voorstelling van een functie op een eindig interval door een Fourierreeks. Eender welke basis van L2 ([a, b]) kan worden gebruikt om een stochastisch proces voor te stellen als een expansie in deze vorm, zoals voor het eerst bestudeerd door Kosambi [25]. Het belang van de (Kosambi–)Karhunen–Loéve-expansie, genoemd naar de Fin Kari Karhunen en de Amerikaan Michel Loève, bestaat er in dat deze de beste basis levert die de mean square error (MSE) op de eindige voorstelling minimaliseert. De coëfficiënten ξn in deze spectrale voorstelling zijn ongecorreleerde toevalsgetallen en de orthogonale basisfuncties fn worden afgeleid uit de covariantiefunctie (2.3). Beschouwen we het stochastisch proces Z(x; ω), x ∈ D, ω ∈ Ω met verwachte ¯ waarde over alle mogelijke realisaties Z(x) en covariantiefunctie C(x, y) (dezelfde als in (2.3)). De covariantiefunctie, ook wel kernel genoemd, is een continue symmetrische positief semidefiniete functie die de covariantie van twee waarden van het stochastische veld Z op posities x en y uitdrukt: C(x, y) , cov(Z(x; ω), Z(y; ω)). Met deze covariantiefunctie is een integraaltransformatie T geassocieerd van de volgende vorm: (T φ)(y) = 1

Z x2

C(x, y)φ(x)dx

x1

Merk op dat de KL-expansie kan worden gezien als aan speciaal geval van polynomial chaos expansion met maximale graad van de veelterm K = 1 [10].

16

2.2. Samplen uit een stochastisch veld waarbij C de kern (kernel) van de transformatie wordt genoemd en φ ∈ L2 [x1 , x2 ]. Omdat T een lineaire operator is, is het zinvol om te spreken over diens eigenwaarden en eigenfuncties. Volgens het Theorema van Mercer bestaat er een orthonormale basis {fn }n voor L2 [x1 , x2 ]. De functies fn zijn de eigenfuncties van de lineaire operator T horende bij de set van eigenwaarden {θn }n van operator T . Merk op dat deze eigenwaarden θn niet negatief zijn. De eigenwaarden θn en eigenfuncties fn voldoen aan het eigenwaardenprobleem Z D

C(x, y)fn (x)dx = θn fn (y).

(2.5)

Verder geldt dat er een spectrale decompositie bestaat van de vorm C(x, y) = P∞ n=0 θn fn (x)fn (y). Het zijn precies deze eigenwaarden en eigenfuncties die gebruikt worden in de expansie (2.4). Men kan aantonen dat de toevalsgetallen ξn (ω) in de KL-expansie (2.4) moeten voldoen aan [16] E[ξn (ω)] = 0 en

E[ξn (ω)ξm (ω] = δnm .

Enkele eigenschappen van de KL-expansie zijn [16]: • Het coördinatensysteem gedefinieerd door de eigenfuncties van de kernel is optimaal in de zin dat het de MSE op de eindige voorstelling van het stochastisch proces Z minimaliseert (in de L2 -norm). • Wanneer het proces Z een Gaussiaans proces is met gegeven covariantiefunctie C, dan zijn de getallen ξn (ω) in de KL-expansie normaal verdeeld. In het geval p = 1 in (2.3) en D = [0, 1] kan men analytische uitdrukkingen voor de eigenwaarden θn en eigenfuncties fn vinden: fn (x) = An (sin(wn x) + λwn cos(wn x))

n = 1, 2, 3...

.

(2.6)

Hierin is wn de n-de van nul verschillende positieve oplossing van de transcendente vergelijking tan w = λ22λw . De constante An is een normalisatieconstante zodat w2 −1 kfn kL2 [0,1] = 1. De overeenkomstige eigenwaarden zijn θn =

2λσ 2 . λ2 wn2 + 1

(2.7)

De eigenwaarden staan weergegeven in Figuur 2.3. In die figuur wordt ook het belang van de correlatielengte λ in (2.3) aangetoond. Zie Appendix A voor de afleiding van deze uitdrukkingen. Voor p = 2 en andere vormen van de covariantiefunctie moet men meer geavanceerde technieken zoals circulant embedding gebruiken. Voor de eigenschappen van stochastische velden in twee dimensies, kan men gebruik maken van een tensorproduct: θn2D = θi1D θj1D , n n

fn2D (x) = fi1D (x1 )fj1D (x2 ) voor in , jn ∈ N. n n

(2.8) 17

2. Multilevel Monte-Carlo voor elliptische SPDEs λ = 0.01 10−1

λ = 0.1 λ=1

θn

10−3

10−5

10−7 100

101

102

103

n Figuur 2.3: Eigenwaarden θn in één dimensie voor λ = 0.01, 0.1, 1 en σ 2 = 1 ten opzichte van de index n. In de praktijk moet men de oneindige som (2.4) afbreken na een eindig aantal termen mKL : Z ≈ ZmKL . Een benadering van de covariantiefunctie C(x, y) opgesteld met mKL = 4 termen staat weergegeven in Figuur 2.2. Om te bepalen hoeveel termen vereist zijn om een zekere nauwkeurigheid te bereiken, kan men volgende eigenschappen gebruiken: • de grootte van de eigenwaarde θn neemt kwadratisch af met n: θn h O(n−2 );

• wanneer de correlatielengte λ < diam(D) is er een “plateau” waar de eigenwaarden nog geen asymptotisch gedrag vertonen: alle eigenfuncties horende bij deze eigenwaarden dragen evenveel bij tot de voorstelling van het stochastische veld (zie Figuur 2.3);

• en, tenslotte,

∞ X

θn = σ

n=1

2

Z

(2.9)

dx.

D

Omdat de getallen ξn (ω) ongecorreleerd zijn volgt hieruit dat 

E kZ − ZmKL k2 = E k h

i

∞ X n=mKL +1

 p

θn fn ξn (ω)k2  =

∞ X n=mKL +1

θn

Omwille van de eerste observatie dat θn h O(n−2 ) volgt dat de RMSE door het −1/2 afbreken van de KL-expansie O(mKL ) is. Een zeer groot aantal termen mKL is dus vereist om een aanvaardbare relatieve fout te bereiken.

2.3

Eindige-volumemethode

In §2.2 hebben we een manier opgesteld om een sample te nemen van de hydraulische conductiviteitstensor k(x; ω). Voor een gegeven sample moet nu de oplossing van een 18

2.3. Eindige-volumemethode deterministische PDE worden berekend. In de context van de numerieke stromingsleer (computational fluid dynamics of CFD) wordt vaak een eindige-volumemethode gebruikt [15]. De behoudswetten worden dan toegepast op een controlevolume (CV) zodat dat de numerieke flux wordt behouden [35]. Het verschil met de eindigedifferentiemethodes is dat het schema wordt uitgeschreven voor de flux in plaats van voor de operator zelf. In wat volgt bespreken we de eindige-volumemethode in één en twee dimensies. Het is belangrijk om op te merken dat de gebruikte discretisatietechniek geen rol speelt in het MLMC-algoritme; men kan in principe eender welke methode gebruiken voor het oplossen van een deterministische PDE.

2.3.1

Eindige-volumemethode in één dimensie

De PDE (2.2) herleidt zich tot een ODE van de vorm d d k(x; ω) p(x) = f (x) − dx dx

x ∈ [0, 1].

(2.10)

Beschouw dan de volgende discretisatie: verdeel [0, 1] in m even grote deelgebieden Di van lengte 1/m met i = 1...m. Die deelgebieden noemt men ook wel cellen of controlevolumes (CV’s). Het centrum van elk deelgebied is een punt, dat we noteren als xi . Integreren van de PDE over elke cel geeft gelijkheid Z Di

−∇ · (k∇p) =

Z

f

Di

1 ≤ i ≤ m.

De behoudswet (2.2) wordt dus uitgeschreven voor elk CV Di afzonderlijk. Door het toepassen van de divergentiestelling van Gauss kan de integraal over het volume Di als een oppervlak-integraal over de rand van het gebied worden geschreven. Elk van de bijdragen op de rand kunnen dan één voor één worden benaderd. Z ∂Di

−k∇p · n =

Z

f

Di

1 ≤ i ≤ m.

(2.11)

Merk op dat de integraal over de rand van het gebied in het ééndimensionale geval zich herleidt tot een evaluatie van de integrand in de eindpunten van het interval. Veronderstel ki en fi de waarde van k en f in het centrum xi van het deelgebied Di . Dan zoeken we een benadering pi voor de waarde van de druk p in xi . Het rechterlid van (2.11) kan worden benaderd via de middelpuntsregel als Z x 1 i+ 2

xi− 1

f (x)dx = (xi+ 1 − xi− 1 )fi ≈ 2

2

fi . m

(2.12)

2

Het linkerlid is de som van de bijdrage op de rand tussen Di−1 en Di , en Di en Di+1 . Om k op de rand te benaderen gebruikt men vaak het harmonisch gemiddelde k¯i+ 1 = 2

2 1 1 + ki ki+1

=

2ki ki+1 . ki + ki+1 19

2. Multilevel Monte-Carlo voor elliptische SPDEs Voor het benaderen van ∇p · n, de naar buiten gerichte normale afgeleide, op de rand kan men een centrale differentie gebruiken: pi − pi−1 , xi − xi−1

∇p · n|∂i−1,i = −

en

∇p · n|∂i,i+1 =

De gediscretiseerde ODE (2.10) wordt dan − −k¯i− 1

2

pi − pi−1 ∆x

+ k¯i+ 1

2

pi+1 − pi ∆x

pi+1 − pi . xi+1 − xi =

fi , m

en, met ∆x = 1/m, −k¯i− 1 pi−1 + Σi pi − k¯i+ 1 pi+1 = 2

2

fi , m2

waarin Σi = (k¯i− 1 + k¯i+ 1 ). In de numerieke experimenten nemen we voor het 2 2 ééndimensionale geval randvoorwaarden van het Dirichlet-type aan: p(x = 0) = 1 en p(x = 1) = 0. Voor de linkerrand vervangen we het harmonisch gemiddelde door de waarde van k in de eerst cel en benaderen we de fluxterm door −(p1 − p(x = 0))/(∆x/2). Een analoge redenering geldt voor de linkerrand. Uitgeschreven in matrixvorm moet men dus voor elk sample het volgende symmetrische tridiagonale stelsel oplossen: Σ1  ¯ −k          

2− 21

−k¯1+ 1 2 Σ2 −k¯3− 1



−k¯2+ 1 2 Σ3

2

−k¯3+ 1 2 .. . −k¯

m−1− 21

Σm−1 −k¯ 1 m− 2

  p1   p   2      p3   .  =  .   .    pm−1  −k¯m−1+ 1   2

Σm

pm

f1 /m2 + 2k(x1 )p(x = 0)   f2 /m2     2 f3 /m     ..     . 



(2.13)

fm−1 /m2 fm /m2 − 2k(xm )p(x = 1)

 

 

Voor dergelijke stelsels kan men het zeer efficiënte Thomasalgoritme gebruiken met een tijdscomplexiteit O(m) [35].

2.3.2

Eindige-volumemethode in twee dimensies

Voor het tweedimensionale geval moet men PDE (2.2) oplossen. Beschouw daartoe de discretisatie van D = [0, 1]2 in m2 gelijke vierkante cellen Di,j = 20

"

i−1 i , × m m !

j−1 j , m m

!#

, met i, j = 1 . . . m.

2.3. Eindige-volumemethode Het centrum van elk deelgebied noteren we als xi,j . De lengte van elke cel is ∆x en de oppervlakte ∆x2 . Opnieuw integreren we (2.2) over elke cel Di,j , Z Di,j

−∇ · (k∇p) =

Z

f

1 ≤ i, j ≤ m.

f

1 ≤ i, j ≤ m.

Di,j

Toepassen van de divergentiestelling levert Z ∂Di,j

−k∇p · n =

Z Di,j

(2.14)

We noteren met ki,j en fi,j de waarde van k en f in het centrum xi,j van het deelgebied Di,j . Dan zoeken we een benadering pi,j voor de waarde van de druk p in het centrum xi,j van elke cel. Het rechterlid van (2.14) kan worden benaderd f als fi,j ∆x2 , of nog: mi,j2 . Het linkerlid is de som van de bijdragen op de linker- en rechterkant en boven- en onderkant van elk deelgebied Di,j . Opnieuw gebruiken we het harmonisch gemiddelde k¯ om de doorlatendheidstensor op de rand te benaderen. Voor de rand tussen Di,j en Di+1,j , wordt dit 2ki,j ki+1,j k¯i+ 1 ,j = 2 ki,j + ki+1,j en analoog voor de andere drie bijdragen. Voor het benaderen van de naar buiten gerichte normale afgeleide ∇p · n op de rand kan men opnieuw een centrale differentie gebruiken. Voor de rand tussen xi,j en xi+1,j wordt dit bijvoorbeeld ∇p · n|∂

i+ 1 2 ,j

=

pi+1,j − pi,j . xi+1,j − xi,j

Uiteindelijk bekomt men een analoge discretisatie als (2.3.1) [8] pi+1,j − pi,j pi,j − pi−1,j ∆x − k¯i− 1 ,j ∆x + 2 ∆x ∆x pi,j+1 − pi,j pi,j − pi,j−1 fi,j ∆x − k¯i,j− 1 ∆x = 2 , 2 ∆x ∆x m

− k¯i+ 1 ,j

2

k¯i,j+ 1 2

of, met ∆x = 1/m, fi,j − k¯i− 1 ,j pi−1,j − k¯i+ 1 ,j pi+1,j + Σi,j pi,j − k¯i,j− 1 ,j pi,j−1 − k¯i,j+ 1 pi,j+1 = 2 . (2.15) 2 2 2 2 m Het symbool Σi,j staat nu voor k¯i− 1 ,j + k¯i+ 1 ,j + k¯i,j− 1 + k¯i,j+ 1 ). Voor de berekening 2 2 2 2 van de waarde pi,j in cel Di,j worden de vier naburige punten gebruikt. De bijdragen kunnen worden geordend als een molecule of stencil [35] van de vorm

. 21


Figuur 2.4: De structuur van de coëfficiëntenmatrix, corresponderende met de eindige volumediscretisatie, vergelijking (2.15).

In de numerieke experimenten nemen we voor het tweedimensionale geval de klassieke randvoorwaarden aan van een zogenaamde flow cel [21]: Dirichletrandvoorwaarden links en rechts en randvoorwaarden van het Neumanntype boven- en onderaan: p(x1 = 0) = 1,

p(x1 = 1) = 0,

∂p = 0, ∂n x2 =0

∂p = 0. ∂n x2 =1

Een voorgeschreven flux kan men gemakkelijk aanleggen door de bijdrage op de rand te vervangen door de waarde van de fluxterm vermenigvuldigd met ∆x. De behandeling van de Dirichletrandvoorwaarden gebeurt analoog als bij het ééndimensionale geval. De structuur van het resulterende lineaire systeem staat weergegeven in Figuur 2.4. Het linkerlid in (2.15) geeft aanleiding tot een ijle matrix met een zeer slecht conditiegetal. Het kan efficiënt worden opgelost via een directe methode of met behulp van iteratieve methodes. Verderop in deze tekst gebruiken we een ijle Choleskyfactorisatie A = LLT , gevolgd door het oplossen van twee stelsels met een driehoeksmatrix als coëfficiëntenmatrix.

2.4

Resultaten

In wat volgt passen we de MLMC-methode uit Hoofdstuk 1 toe op het modelprobleem uit §2.1. De sequentie van benaderingen (1.5) kiezen we als M` = m` = (m0 2` )d met m0 de roosterafstand op het ruwste rooster en d = {1, 2} de dimensie van het probleem. De roosterafstand op level ` is dan h` = m−1 ` , waarbij de roosterafstand bij elke verfijning gehalveerd wordt. In eerste instantie beschouwen we het ééndimensionale geval d = 1. We maken een gefundeerde analyse van de winst aan rekenkost ten opzichten van standaard Monte-Carlo-technieken. Tot slot bestuderen we ook het tweedimensionale geval. 22

2.4. Resultaten

2.4.1

Resultaten voor het ééndimensionale geval

De ééndimensionale probleemstelling is eerder al geschetst; PDE (2.2) herleidt zich tot ODE (2.10) en de randvoorwaarden kiezen we als p(0) = 1 en p(1) = 0. De eindfunctie, waarvan we de stochastische karakteristieken wensen te berekenen, is de uitstroming op x = 1, ∂p . (2.16) Q = −k ∂x x=1 We berekenen deze uitstroom door middel van een eerste-orde differentie op de rechterrand: Q = h2` k(xm )p(xm ). De covariantiefunctie C(x, y) van de doorlatendheidstensor k(x; ω) kiezen we van de vorm (2.3) met variantie σ 2 = 1 en correlatielengte λ = 0.3. Het aantal termen in de KL-expansie, mKL , is 800 en de ruwste roosterafstand m0 = 16. Figuur 2.5 bovenaan links toont de variantie van de eindfunctie Q` op elk level ` en het verschil Y` = Q` − Q`−1 gebaseerd op N =1e5 samples. Merk op dat V[Q` ] inderdaad constant is, zoals verwacht, terwijl V[Y` ] een sterk dalende functie is. De helling van de lijn V[Y` ] is ongeveer -2, waaruit volgt dat V[Y` ] . M`−β en dus β ≈ 2. Analoog, en kijken naar Figuur 2.5 bovenaan rechts, de helling van de lijn E[Y` ] is ongeveer -1.75, waaruit volgt dat E[Y` ] . M`−α en dus dat α ≈ 1.75. Figuur 2.6(a) toont de reële uitvoeringstijd in seconden voor het nemen van N = 1000 samples in één dimensie op verschillende levels `. Dit bevestigt dat in het ééndimensionale geval γ ≈ 1. De voorwaarden (i) − (iii) van Theorema 1 zijn dus voldaan met β > γ. Dat wil zeggen dat er voor eender welke tolerantie op de RMSE, hoe klein ook, een bovengrens van O(2 ) is op de kost van het MLMC algoritme. Dit wordt bevestigd in Figuur 2.5, onderaan rechts. Daar wordt het product van 2 met de rekenkost uitgezet als functie van . Zoals voorspeld door de theorie, krijgen we een (bijna) constante waarde. Figuur 2.5 onderaan links toont het aantal samples op elk level ` voor verschillende toleranties op de RMSE in het MLMC-algoritme. Merk op dat de grotere variantie op level 0 wordt weerspiegeld in het grote aantal samples benodigd op dat level. Voor de rekencomplexiteit is dit een goede zaak. Level 0 stelt immers het ruwste rooster voor, waarop samples van PDE (2.2) het goedkoopst zijn. Figuur 2.5 onderaan rechts toont een vergelijking van de kost van Monte-Carlo-simulatie en multilevel Monte-Carlo-simulatie. Deze standaardkosten werden berekend als C MLMC = N0 M0 +

L X `=1

N`

M` + M`−1 M0

(2.17)

voor multilevel Monte-Carlo en C MC =

L X

N`∗ M`

(2.18)

`=0

voor gewone Monte-Carlo waarin N`∗ = 2−2 V[Q` ] (zie [17]). Op die manier wordt hetzelfde heuristische stopcriterium op elk level gebruikt. De kost voor MLMC is 23

2. Multilevel Monte-Carlo voor elliptische SPDEs 0

−5

−5

log2 V[ · ]

log2 |E[ · ]|

0

−10

−10

−15

−20

−15

Q` Q` − Q`−1 0

2

−20

4

Q` Q` − Q`−1 0

2

` 1010

= = = =

103

0.001 0.0005 0.0001 0.0005

2 cost

N`

108

4

`

106

MC MLMC

102

104 0

1

2

`

3

4

101

10−4

10−3

Figuur 2.5: Performantie van MLMC-simulatie voor het modelprobleem (2.1) met λ = 0.3, σ 2 = 1, s = 2, mKL = 800 en m0 = 16 in één dimensie. De eindfunctie is de ∂p uitstroom −k ∂x op x = 1. niet alleen veel lager dan voor MC, ze neemt ook nauwelijks toe met afnemende tolerantie . In wat volgt bestuderen we de winst aan rekenkost van het MLMC-algoritme meer gedetailleerd. We nemen de ruimtelijke discretisatie op het fijnste rooster constant, ML = 16 · 24 = 256, dus L = 4. In Figuur 2.7(a) analyseren we hoe de kost van de methode toeneemt naarmate de tolerantie δ op de standaarddeviatie van de schatter daalt. De standaarddeviatie van de schatter is de vierkantswortel uit de stochastische fout (de variantie van de schatter) in (1.9). Een√tolerantie δ op de standaarddeviatie komt dus overeen met een tolerantie = 2δ op de RMSE (1.3). De geschatte ruimtelijke discretisatiefout op het fijnste rooster is 1 ∼ E[QM − Q] = E[QL − QL−1 ] M α −1 ≈ 1.0616e-04, zie (1.17). De eindfunctie is

∂p opnieuw de uitstroom op x = 1, −k ∂x . De winst aan rekenkost door het toevoegen van een multilevelstructuur is duidelijk zichtbaar. De rekenkost van de standaard Monte-

24

2.4. Resultaten 30

CPU-tijd [s]

CPU-tijd [s]

20

15

10

5

0

2

4

ML (a) 1D

6

20

10

tijd [s]

tijd [s]

O(ML )

γ O(ML )

8 ·103

0

0

2

4

ML

6 ·104

(b) 2D

Figuur 2.6: (a) CPU-tijd voor L =0:9, m0 = 16, mKL = 800, λ = σ 2 = 1 en N = 1000 oplossingen van de elliptische PDE in één dimensie met lognormaal verdeelde k uitgezet ten opzichte van het aantal roosterpunten in de discretisatie, ML = mL = m0 2L . (b) CPU-tijd voor L =0:6, m0 = 4, mKL = 1400, λ = σ 2 = 1 en N = 50 oplossingen van de elliptische PDE in twee dimensies met random gesamplede k ten opzichte van het aantal roosterpunten in de discretisatie, ML = (mL )2 = (m0 2L )2 . De sparse direct solver van Matlab voor het oplossen van het stelsel in 2D is CHOLMOD. In dit geval is γ = 1.1436. Carlo-schatter met een tolerantie op de standaardafwijking van 2e-4 is even groot als de rekenkost van de 4-level-schatter met een tolerantie op de standaardafwijking van 1e-4. Dit ligt onder de geschatte discretisatiefout op het fijnste rooster. In Figuur 2.7(b) veronderstellen we de gevraagde tolerantie op de standaardafwijking van de schatter constant, δ =1e-3, en onderzoeken we hoe de kost van de MLMC-methode met een vast aantal levels L =1:4 toeneemt met ML , het aantal roosterpunten in één dimensie. Hier blijkt duidelijk de verbetering van MLMC ten opzichte van standaard MC. De kost voor het schatten van E[Q` ] op een rooster met ML = 32 met standaard Monte Carlo is even groot als de kost voor het schatten van E[Q` ] op een rooster met ML = 128 met de 4-level MLMC-methode. De theoretische winst in rekenkost uit zich ook in de reële uitvoeringstijd, zoals getoond in Figuur 2.8. Deze figuur toont dezelfde resultaten als Figuur 2.7(b), maar met de reële uitvoeringstijd op de y-as, in plaats van de standaardkost.

2.4.2

Resultaten voor het tweedimensionale geval

Tot slot bespreken we ook enkele resultaten van een simulatie in twee dimensies. Net als in het ééndimensionale geval kiezen we de variantie σ 2 = 1 en de correlatielengte λ = 0.3. Het aantal termen in de KL-expansie is mKL = 1400 waarbij de eigenfuncties worden geordend in volgorde van dalende eigenwaarden θn . De ruwste roosterafstand 25


108

10

MC 2-level MC 3-level MC 4-level MC

standaardkost

δ

10−2

−3

107

MC 2-level 3-level 4-level 5-level

MC MC MC MC

106 10−4 104

105

106

107

108

101

109

102

standaardkost

ML

(a)

(b)

CPU tijd [s]

Figuur 2.7: (a) Gevraagde tolerantie op de standaardafwijking δ van de schatter ˆ MLMC , L = 1:4. versus gestandaardiseerde rekenkost voor verschillende schatters Q M,{N` } De horizontale streeplijn is de geschatte discretisatiefout op het rooster. De eindfunc∂p tie is de uitstroom op x = 1, −k ∂x . (b) Gestandaardiseerde rekenkost ten opzichte van het aantal roosterpunten ML voor een vaste tolerantie δ =1e-3 op de maximale standaardafwijking van de MLMC–schatter met een vast aantal levels L =1:5.

MC 2-level 3-level 4-level 5-level

MC MC MC MC

102

101

102

M LL Figuur 2.8: Reële uitvoeringstijd (CPU tijd in seconden) ten opzichte van het aantal roosterpunten ML voor een vaste tolerantie δ =1e-3 op de maximale standaardafwijking van de MLMC–schatter met een vast aantal levels L =1:5.

26

2.5. Besluit is m0 = 8. Als eindfunctie kiezen we de effectieve conductiviteit [21] Q = keff , −

Z 1 0

∂p dx2 , ∂x1 x1 =1

k

(2.19)

de totale uitstroom op x1 = 1. De afgeleide wordt opnieuw benaderd met een eindige differentie en de integraal wordt benaderd door een kwadratuurregel. Omdat de druk p slechts in discrete punten gekend is, kiezen we voor de trapeziumregel [23]. Voor het tweedimensionale geval geldt dat M` = (m` )2 = (m0 2` )d . Men moet voor elk sample een bandstelsel (2.15) oplossen. Figuur 2.9 bovenaan toont de variantie van de eindfunctie Q` op elk level ` en het verschil Y` = Q` − Q`−1 gebaseerd op N =1e6 samples. De helling van de lijn V[Y` ] is ongeveer −3, waaruit volgt dat V[Y` ] . M`β en dus β ≈ 3. Analoog, de helling van de lijn E[Y` ] is ongeveer −0.75, waaruit volgt dat E[Y` ] . M`α en dus dat α ≈ 0.75. Figuur 2.6(b) toont de reële uitvoeringstijd in seconden voor het nemen van N = 50 samples in twee dimensie op verschillende levels `. In het tweedimensionale geval is γ ≈ 1.1436. Figuur 2.9 onderaan links toont opnieuw het aantal samples op elk level ` voor verschillende toleranties op de RMSE in het MLMC-algoritme. Figuur 2.9 onderaan rechts toont een vergelijking van de kost van Monte-Carlo-simulatie en multilevel Monte-Carlo-simulatie. Bemerk opnieuw de verbetering van MLMC ten opzichte van gewone MC, in overeenstemming met Theorema 1.

2.5

Besluit

Voor het modelprobleem van een elliptische partiële differentiaalvergelijking met onzekerheid in de diffusiecoëfficiënt levert het gebruik van de multilevel Monte-Carloschatter een aanzienlijke winst aan rekenkost op ten opzichte van de standaard Monte-Carlo-schatter. Ook voor stochastische gewone differentiaalvergelijkingen (SDEs) is het gebruik van een mulilevelschatter voordelig [17]. Voor PDEs is het voordeel echter nog heel wat groter. Voor een PDE toepassing is het nemen van een sample op een ruw rooster immers vele malen goedkoper dan op een fijn rooster. In dit hoofdstuk toonden we dat voor het samplen uit een stochastisch veld k(x; ω) in het geval p = 1 in (2.3) er analytische uitdrukkingen voor eigenwaarden en eigenfuncties bestaan. De resulterende deterministische PDE kan dan worden opgelost met bijvoorbeeld een eindige-volumemethode. Enkele beknopte numerieke experimenten toonden de effectiviteit en de effciëntie van het algoritme aan, zowel in één als twee dimensies. In de recente onderzoeksliteratuur gaat men op zoek naar nog snellere varianten van het multilevelalgoritme. Door gebruik te maken van variantiereductietechnieken probeert men de variantie van de termen in de telescopische som (1.6) nog verder te verlagen. Enkele voorbeelden zijn MLMC met importance sampling en Continuation Multilevel Monte Carlo (CMLMC), waar de MSE (1.9) bestaat uit de gewogen som van de stochastische fout en de discretisatiefout. In het vervolg van dit eindwerk bestuderen we nog twee technieken om het MLMCalgoritme te versnellen. In hoofdstuk 3 onderzoeken we hoe het algoritme kan worden 27

2. Multilevel Monte-Carlo voor elliptische SPDEs 0

0

−2

log2 V[ · ]

log2 |E[ · ]|

−5

−10

−15

−20

1

2

−6 −8

Q` Q` − Q`−1 0

−4

3

4

−10

5

Q` Q` − Q`−1 0

`

2

3

4

5

`

105 10

1

4

= = = =

105

0.1 0.05 0.02 0.01

MC MLMC

104

N`

2 cost

103 102

103

101 100

0

1

2

3

`

4

5

102 10−2

10−1

Figuur 2.9: Performantie van MLMC-simulatie voor het modelprobleem λ = 0.3, σ 2 = 1, mKL = 1400 en ruwste roosterafstand m0 = 8 in twee dimensies. De eindfunctie is keff (2.19). versneld door parallelle uitvoering en bestuderen we zogenaamde levelafhankelijke schatters. In hoofdstuk 4 gebruiken we quasi-Monte-Carlo-punten voor het samplen uit het stochastisch veld en komen zo tot een multilevel quasi-Monte-Carlo-algoritme (MLQMC).

28

Hoofdstuk 3

Analyse van de MLMC-methode voor het modelprobleem In dit hoofdstuk wordt het modelprobleem, een elliptische partiële differentiaalvergelijking met onzekerheid in de diffusiecoëfficiënt, in meer detail besproken. In een eerste deel bekijken we hoe de MLMC-benadering kan worden geparallelliseerd. Het zal blijken dat, net als gewone Monte-Carlo, de multilevelbenadering efficiënt in parallel kan worden uitgevoerd. Daarna wordt een parameteranalyse uitgevoerd: hoe verandert de oplossing van het modelprobleem in één en twee dimensies als de parameters in de stochastische coëfficiënt worden gewijzigd. Tot slot bespreken we de toepassing van zogenaamde levelafhankelijke schatters; deze laten toe om ook voor heel kleine correlatielengten λ de multilevelstructuur te gebruiken.

3.1

Parallellisatie

Net als gewone Monte-Carlo-simulatie leent multilevel Monte-Carlo zich uitstekend tot parallelle verwerking. Gegeven het aantal samples op elk level kan men elk van deze deterministische problemen onafhankelijk van elkaar oplossen. Er is immers geen communicatie tussen de verschillende processoren vereist wanneer de beschikbare hoeveelheid geheugen voldoende is om elk sample op één processor te berekenen. In wat volgt geven we enkele definities die vaak worden gehanteerd binnen het domein van de parallel computing. Deze laten dan toe om de performantie van het parallelle algoritme te vergelijken met de sequentiële implementatie.

3.1.1

Definities

De meest natuurlijke grootheid om performantie van een parallel algoritme te beschrijven is de speedup. Deze is gedefinieerd als volgt [2]: Definitie 3.1. (speedup) Gegeven Tseq , de sequentiële uitvoeringstijd van een algoritme, en Tpar , de parallelle uitvoeringstijd gebruikmakende van p processoren, dan 29

3. Analyse van de MLMC-methode voor het modelprobleem is de speedup S gegeven door S=

Tseq . Tpar

(3.1)

Ideale speedup betekent dat S = p, immers, dan is de parallelle uitvoeringstijd Tpar gelijk aan de sequentiële uitvoeringstijd Tseq gedeeld door het aantal processoren. Wanneer S < 1 spreekt men van slowdown. Door gebruik te maken van caches is het soms ook mogelijk om S > p te bereiken. Niet alle algoritmes zijn even goed parallelliseerbaar. Om een realistische vergelijk te maken dient men voor Tseq het beste sequentiële algoritme te gebruiken in plaats van de uitvoeringstijd op p = 1 processoren. Op die manier vergelijkt men de beste sequentiële implementatie met de beste parallelle implementatie, die eventueel een ander algoritme gebruikt. De speedup S hangt af van zowel het aantal processoren p als de probleemgrootte n. Men kan het geparallelliseerde algoritme dus vergelijken met de sequentiële implementatie voor zowel een variërend aantal processoren p als een variërende probleemgrootte. Definitie 3.2. (sterke schaalbaarheid) Een algoritme heeft een sterke schaalbaarheid wanneer Tseq S(p) = = O(p). (3.2) Tpar

Definitie 3.3. (zwakke schaalbaarheid) Een algoritme heeft een zwakke schaalbaarheid wanneer Tseq (n) S(n) = = O(1). (3.3) Tpar (n)

In een realistische toepassing is het niet redelijk om sterke schaalbaarheid te verwachten. Men zal immers nooit hetzelfde probleem oplossen met een verschillend aantal processoren. Bovendien is een algoritme niet oneindig parallelliseerbaar: wanneer er evenveel processoren als elementaire taken1 zijn heeft het geen zin om extra processoren te gebruiken. In de praktijk is het aantal processoren een constante, en is de probleemgrootte een functie van dit aantal processoren. Het is daarom wel logisch om zwakke schaalbaarheid van een algoritme te verwachten wanneer n → ∞. 1 een elementaire taak of chunk is een deelprobleem dat klein genoeg is om op meerdere processoren te worden uitgevoerd maar groot genoeg om de overhead van de parallelle uitvoering te overwinnen. Dit noemt men ook wel de granulariteit. Hoe kleiner de granulariteit, hoe groter de potentiële winst door parallellisatie en dus de speedup, maar ook hoe groter de synchronisatie- en communicatiekost.

30

3.1. Parallellisatie

3.1.2

Methode

Men kan het MLMC-algoritme op drie verschillende manieren parallelliseren. 1. Parallellisatie van de PDE-solver: in elk deterministisch sample kan men het ijle stelsel dat moet worden opgelost oplossen met een iteratieve methode. De basis van elke iteratieve methode is de matrix-vector vermenigvuldiging. Door een gepaste matrixpartitionering is het mogelijk om deze vermenigvuldiging in parallel uit te voeren. 2. Homogene parallellisatie: op elk level ` kan men het optimale aantal samples N` in parallel nemen. Deze methode is het eenvoudigst, en heeft als bijkomend voordeel dat elke deterministische oplossing van de PDE op een gegeven rooster in principe even veel tijd in beslag neemt. Op die manier heeft men automatisch een gepaste werkverdeling (load balance). 3. Heterogene parallellisatie: men neemt samples in parallel over de verschillende levels heen. In alle volgende numerieke experimenten opteren we voor de tweede benadering. Een mogelijke beperking voor deze methode is geheugengebruik. Elke processor moet over voldoende geheugen beschikken om één sample op het fijnste rooster onafhankelijk te kunnen berekenen. In de praktijk blijkt dit altijd het geval. Het fijnste rooster L is eerder beperkt door het geheugengebruik voor de opslag van alle eigenfuncties (tot 1400 in het tweedimensionale geval) dan voor het oplossen van de PDE. Voor complexe geometrieën of fijnere roosters kan men gebruik maken van domeindecompositie om deze hindernis te overwinnen. We gebruiken een Matlab-implementatie van het MLMC-algoritme op de computercluster van de KU Leuven. Elke compute node bestaat uit twee 10-core Ivy Bridge 2.8GHz Xeon E5-2680v2 processoren met 64 GByte RAM en 25 MByte level 3 cache.

3.1.3

Performantie

We analyseren de performantie van de parallelle implementatie aan de hand van Figuur 3.1. De linkse grafiek toont de speedup (3.1) ten opzichte van het aantal processoren voor een vaste probleemgrootte (tolerantie) van = 1e-3. Voor het ééndimensionale geval is er een bijna perfecte speedup tot op p = 20 processoren. Voor d = 2 presteert het algoritme minder optimaal. Het MLMC-algoritme is dus bijna sterk schaalbaar, en dit is ook wat men kan verwachten: er is geen overhead door synchronisatie of communicatie omdat alle samples onafhankelijk van elkaar kunnen worden berekend. De rechtse grafiek toont speedup ten opzichte van de probleemgrootte. Hier speelt de tolerantie op de RMSE de rol van probleemgrootte, het is immers de enige instelbare parameter voor de hoeveelheid werk in het MLMC-algoritme. Merk op hoe zowel voor d = 1 als d = 2 er geen perfecte speedup wordt bereikt. Wanneer de 31

3. Analyse van de MLMC-methode voor het modelprobleem 20

25 d=1

d=1

d=2

S()

S(p)

10

15

10

5

0

d=2

20

15

5

0

5

10

15

0

20

10−2

p

10−3

10−4

(a)

(b)

Figuur 3.1: (a) Speedup versus aantal processoren. (b) Speedup versus probleemgrootte, hier . speedup zich zou stabiliseren voor nog kleinere toleranties dan is het MLMC-algoritme zwak schaalbaar.

3.2

Parameteranalyse van het modelprobleem

De twee grootheden die het stochastische veld k(x; ω) karakteriseren zijn de correlatielengte λ en variantie σ 2 . Afhankelijk van het type poreus materiaal in het modelprobleem kan men een verschillend stochastisch model met een verschillende correlatielengte en variantie fitten voor k. Men kan zich nu afvragen wat het effect is van zo’n verandering op de eindfunctie. Hieronder beschouwen we opnieuw het modelprobleem in één en twee dimensies, maar met verschillende waarden voor de parameters in het stochastisch model.

3.2.1

Het 1D-probleem

In het ééndimensionale geval zijn we geïnteresseerd in de uitstroom Q = −k

∂p . ∂x1 x1 =1

(3.4)

We beschouwen daartoe volgende combinaties van correlatielengte λ en variantie σ 2 van het stochastisch veld k(x; ω): A1 B1 C1 D1 E1 32

λ= λ= λ= λ= λ=

0.30, 3.00, 0.03, 0.30, 0.30,

σ2 σ2 σ2 σ2 σ2

= = = = =

1.00 1.00 1.00 0.50 2.00

mKL mKL mKL mKL mKL

= 800 = 250 = 2500 = 500 = 1200

3.2. Parameteranalyse van het modelprobleem Hoe kleiner de correlatielengte en hoe groter de variantie hoe meer termen mKL er nodig zijn in de KL-expansie van stochastische veld en dus hoe moeilijker het probleem. Voor het standaardprobleem λ = 0.3 en σ 2 = 1 is het aantal termen in de KL-expansie zo gekozen dat de kleinste eigenwaarde (2.7) grootte-orde 1e-6 is. Voor de andere gevallen hanteren we een gelijkaardig criterium, het aantal termen is eveneens weergegeven in bovenstaande tabel. Bemerk reeds het effect van een verandering van de variantie σ 2 op het aantal termen dat moet worden meegenomen in de ontwikkeling. De roosterafstand op het ruwste rooster is steeds h0 = 1/16. Figuur 3.2 toont de verwachte waarde en variantie van de eindfunctie Q` en het verschil Y` = Q` − Q`−1 op elk level voor elk van bovenstaande gevallen. Door het variëren van de correlatielengte verandert de verwachte waarde van de eindfunctie. Hoe kleiner de correlatielengte, hoe moeilijker het probleem. Dit is zichtbaar in Figuur 3.2(a). De variantie op level ` = 5 voor λ = 0.03 is ongeveer even groot als de variantie op level ` = 2 voor λ = 3. Er zullen dus veel meer samples nodig zijn op elk level om de variantie van de schatter onder de gevraagde tolerantie 2 /2 te brengen. Een analoge redenering geldt voor de onzuiverheid (bias) in Figuur 3.2(b). Naarmate λ stijgt zijn er steeds meer levels nodig om de discretisatiefout kleiner dan √ / 2 te maken. Door het variëren van de variantie σ 2 verandert de verwachte waarde van de oplossing slechts weinig. Hoe kleiner de variantie van het veld k(x, ω), hoe kleiner de variantie van de oplossing. Naarmate de variantie van het veld stijgt wordt het probleem moeilijker: er zijn meer samples nodig op elk level om de stochastische fout te √ beperken en er zijn meer levels nodig om de discretisatiefout onder de gevraagde / 2 te brengen. Bemerk hoe de variatie van σ 2 een eenvoudige verschuiving van de variantie van Q` en Y` teweeg brengt. De parameters α en β in Theorema 1 zijn dus probleemafhankelijk: ze hangen af van de vergelijking (2.2) en de beschouwde eindfunctie (3.4), maar niet van de gebruikte parameters in het stochastische veld k(x; ω). Met α ≈ 1.75 en β ≈ 2 wil 2 dat zeggen dat we voor elke combinatie van λ en σ dezelfde asymptotische kost ˆ MLMC h O(2 ) mogen verwachten, evenwel met een andere constante. C Q M,{N` } Tabel 3.1 toont de resultaten voor een parameteranalyse van het modelprobleem met d = 1 (één dimensie). Voor elk van de beschouwde gevallen (A1 -E1 ) berekenen we de gemiddelde waarde en variantie van de eindfunctie voor verschillende toleranties op de RMSE. De tabel toont eveneens de variantie van de schatter en de geschatte onzuiverheid (in het kwadraat). Beide termen moeten kleiner zijn dan 2 /2. De laatste kolom toont de effectieve rekentijd in seconden. Alle simulaties werden uitgevoerd met behulp van de parallelle implementatie uit §3.1.2 met p = 20 processoren. Bemerk hoe de variantie van de schatter telkens zéér dicht tegen de ideale 2 /2 ligt. Dit bewijst de effectiviteit van het optimalisatieprobleem (1.16). Door het veranderen van de parameters van het probleem verandert ook het conditiegetal van het stelsel, afgeleid uit een eindige-volumediscretisatie, dat moet worden opgelost in elk deterministisch sample van de elliptische SPDE. Tabel 3.2 toont een vergelijking van de grootte-orde van het conditiegetal van het stelsel voor de verschillende parametersets (A1 -E1 ) in Tabel 3.1. Het conditiegetal stijgt naarmate 33

3. Analyse van de MLMC-methode voor het modelprobleem

variantie

variantie schatter

geschatte onzuiverheid2

tijd [s]

9.495741e-01 5.109685e-01 4.999973e-07 1.578444e-08 3.681375e+01

5e-4

9.500669e-01 5.147396e-01 1.249999e-07 1.514970e-08 1.323067e+02

1e-4

9.502175e-01 5.136264e-01 5.000000e-09 7.474089e-10 3.383141e+03

5e-5

9.502280e-01 5.136750e-01 1.250000e-09 7.967480e-10 1.337198e+04

1e-3

1.490367e+00 3.236540e+00 4.999997e-07 3.340251e-09 1.220919e+02

5e-4

1.489540e+00 3.224270e+00 1.250000e-07 7.057171e-10 4.952275e+02

1e-4

1.489787e+00 3.230410e+00 5.000000e-09 3.262991e-10 9.603357e+03

5e-5

1.489816e+00 3.230254e+00 1.250000e-09 3.284747e-10 1.177745e+04

1e-3

6.787930e-01 9.053706e-02 4.999952e-07 4.712891e-08 1.890904e+01

5e-4

6.787325e-01 9.057728e-02 1.249999e-07 4.254384e-08 6.355610e+01

1e-4

6.791906e-01 9.070403e-02 4.999999e-09 2.998931e-09 1.694623e+03

5e-5

6.792615e-01 9.073624e-02 1.250000e-09 2.057654e-10 7.141812e+03

1e-3

9.703343e-01 2.340763e-01 4.999956e-07 6.368884e-09 1.681062e+01

5e-4

9.702242e-01 2.333767e-01 1.249996e-07 3.361971e-09 5.390952e+01

1e-4

9.701911e-01 2.333994e-01 4.999999e-09 2.173015e-09 1.282075e+03

5e-5

9.700108e-01 2.332978e-01 1.250000e-09 1.387852e-10 5.262923e+03

1e-3

9.332392e-01 1.318862e+00 4.999996e-07 6.993698e-08 5.604143e+01

5e-4

9.324719e-01 1.315551e+00 1.250000e-07 7.693208e-08 2.156297e+02

1e-4

9.321773e-01 1.315847e+00 5.000000e-09 3.693642e-10 5.947461e+03

5e-5

9.321329e-01 1.315319e+00 1.250000e-09 3.569637e-10 2.357207e+04

λ = 3, σ 2 = 1

λ = 0.3, σ 2 = 1

1e-3

λ = 0.3, σ 2 = 0.5 λ = 0.03, σ 2 = 1

gemiddelde waarde

λ = 0.3, σ 2 = 2

Tabel 3.1: Parameteranalyse van het modelprobleem in één dimensie.

34

3.2. Parameteranalyse van het modelprobleem geval

`=0

`=1

`=2

`=3

`=4

`=5

A1 B1 C1 D1 E1

102 102 102 102 102

103 102 103 102 103

103 103 104 103 104

104 104 104 104 104

105 104 105 104 105

105 105 106 105 106

Tabel 3.2: Grootte-orde van het conditiegetal van het stelsel in de eindigevolumemethode in één dimensie. ` toeneemt, en hoe moeilijker het probleem (C1 en E1 ) hoe sneller de stijging.

35


0

0

−10

log2 E[| · |]

log2 V[ · ]

−10

−20

−20

Q` Y` λ=3

−30

Q` Y` λ=3

−30

λ = 0.3

λ = 0.3 λ = 0.03

λ = 0.03 −40

0

1

2

3

4

5

−40

0

1

2

4

5

3

4

5

`

`

log2 E[| · |]

0

log2 V[ · ]

0

−10

−10

Q` Y`

−20

Q` Y`

−20

σ 2 = 0.5

σ 2 = 0.5

σ2 = 1

σ2 = 1

σ2 = 2 −30

3

0

1

σ2 = 2 2

3

`

4

5

−30

0

1

2

`

Figuur 3.2: Parameteranalyse van het modelprobleem met d = 1. Variantie (links) en verwachte waarde (rechts) van de eindfunctie Q` en het verschil Y` = Q` − Q`−1 op elk level voor verschillende λ met σ 2 = 1 (bovenaan) en verschillende σ 2 met λ = 0.3 (onderaan). De resultaten zijn gebaseerd op het nemen van N =1e6 samples op elk level.

36

3.2. Parameteranalyse van het modelprobleem

3.2.2

Afleiden van een kansdichtheidsfunctie

Uit alle voorgaande experimenten blijkt dat de eindfunctie (3.4) net als het stochastisch veld k(x; ω) lognormaal verdeeld is. Een lognormale verdeling heeft een kansdichtheidsfunctie fX van de vorm 1 (log x − θ)2 fX (x) = √ exp − 2φ2 x 2πφ "

#

,

x>0

immers, het logaritme van de lognormale verdeling is normaal verdeeld. De parameters θ en φ zijn de verwachtingswaarde en respectievelijk standaardafwijking van de natuurlijke logaritme van de variabele X. Voor deze verdeling geldt dat E[x] = eθ+ 2 φ 1

2

V[x] = eφ − 1 e2θ+φ

2

= eφ − 1

2

2

E [x]

h

i2

.

Gegeven dus de verwachte waarde en variantie van de eindfunctie dan kan men de parameters θ en φ van de lognormale verdeling schatten door het oplossen van een stelsel van twee niet-lineaire vergelijkingen in twee onbekenden met behulp van bijvoorbeeld de methode van Newton-Raphson. Figuur 3.3 toont de kansdichtheidsfunctie voor het ééndimensionale geval voor elk van de gevallen in bovenstaande tabel. Voor een lognormaal verdeeld veld met parameters λ = 0.3 en σ 2 = 1, bijvoorbeeld, is de eindfunctie (3.4) lognormaal verdeeld met parameters θ ≈ −0.2669 (gemiddelde) en φ2 ≈ 0.4389 (variantie).

3.2.3

Het 2D-probleem

In de literatuur zijn een aantal typische eindfuncties Q voorhanden voor het modelprobleem (2.2) in twee dimensies [38]. We bespreken enkele veel voorkomende gevallen. • Puntevaluatie: men beschouwt de druk p(x∗ ; ω) op een bepaald punt x∗ , vaak gekozen zodanig dat het niet samenvalt met een roosterpunt x∗ ⊂ D = [0, 1]2 . In een ruimtelijke dimensie d > 1 is de druk p vaak onbegrensd in de Sobolevruimte H 1 (D). Men regulariseert dit type functionaal dan ook vaak door het te benaderen door een lokaal gemiddelde p(x∗ ; ω) ≈

1 |D∗ |

Z D∗

p(x; ω)dx

waarin D∗ een deelgebied is van D die x∗ omvat. • Flux: net als een puntevaluatie van de druk kan men ook een puntevaluatie van de flux −k∇p benaderen door een lokaal gemiddelde.

Naast lineaire functionalen kan men ook geïnteresseerd zijn in niet-lineaire functionalen: 37

3. Analyse van de MLMC-methode voor het modelprobleem 2

1.2 λ = 0.3, σ 2 = 1

λ = 0.3, σ 2 = 1

2

λ = 3, σ = 1

λ = 0.3, σ 2 = 2

1

2

λ = 0.3, σ 2 = 0.5

λ = 0.03, σ = 1

1.5

fX (x)

fX (x)

0.8 1

0.6 0.4

0.5 0.2 0

0

1

2

3

4

0

5

0

1

2

x

3

4

5

x

Figuur 3.3: Kansdichtheidsfunctie fX van de eindfunctie voor elk van de testgevallen (A1 -E1 ). De eindfunctie (3.4) is telkens lognormaal verdeeld, waarbij de parameters van de verdeling kunnen bepaald worden door het oplossen van een stelsel van twee niet-lineaire vergelijkingen. • Hogere-orde momenten: de meest voor de hand liggende niet-lineaire functionaal is het hogere-orde momenten van een lineaire functionaal: Mk (p) = E[(p)k ] is het k-de orde moment van p(x; ω). • Norm: men beschouwt de norm van het drukveld, kp(x)kL2 (D) . • Uitstroming: men is geïnteresseerd in de totale uitstroom uit (een deel van) de rand Γ van het gebied D. In (2.19), bijvoorbeeld, beschouwden we de flux doorheen de rechterwand van het gebied. In wat volgt beschouwen we voor het tweedimensionale geval een puntevaluatie in x∗ =

3 1 3 1 + , + 4 512 4 512

(3.5)

in plaats van de uitstroom keff . Dit is een roosterpunten dat, met h0 = 1/8, voor ` ≤ 5 nooit samenvalt met een roosterpunt. We opteren hier voor een andere eindfunctie dan in §2.4.2, omwille van de trage convergentie van de verwachte waarde E[|Y` |] [7]. Net als in §3.2.1 beschouwen we enkele combinaties van λ en σ 2 om te onderzoeken wat het effect is van het model voor het stochastisch veld k(x; ω) op de resultaten: A2 B2 C2 D2 E2 38

λ= λ= λ= λ= λ=

0.30, 0.90, 0.10, 0.30, 0.30,

σ2 σ2 σ2 σ2 σ2

= = = = =

1.00 1.00 1.00 0.50 2.00

mKL mKL mKL mKL mKL

= 1400 = 500 = 2900 = 600 = 2800

3.3. Levelafhankelijke schatters Voor het standaardprobleem λ = 0.3 en σ 2 = 1 in twee dimensies is het aantal termen in de KL-expansie zo gekozen dat de kleinste eigenwaarde grootte-orde 1e-5 is. Merk op dat in twee dimensies θn2D = θi1D · θj1D . Voor de andere gevallen hanteren we een n n gelijkaardig criterium, het aantal termen is eveneens weergegeven in bovenstaande tabel. Een analoge analyse als voor Figuur 3.2 leert dat de verwachte waarde van de eindfunctie (3.5) slechts weinig verandert bij de variantie van de parameters λ en σ 2 . Naarmate de correlatielengte daalt of de variantie stijgt neemt ook de variantie en verwachte waarde van Y` toe. Dat betekent dat men meer samples en/of meer levels moet toevoegen om aan een gegeven tolerantie op de RMSE te voldoen. Het probleem wordt dus moeilijker, iets wat ook zichtbaar is in Tabel 3.3. De parameters α en β in Theorema 1 zijn opnieuw probleemafhankelijk: ze hangen af van de vergelijking en de beschouwde eindfunctie maar niet van de gebruikte parameters in het model voor k(x; ω). Tabel 3.3 toont de resultaten voor een parameteranalyse van het modelprobleem met d = 2 (twee dimensies). Voor een gegeven RSM- berekenen we de gemiddelde waarde en variantie van de druk in x∗ voor elke combinatie A2 t.e.m. E2 . De derde en vierde kolom tonen de variantie van de schatter en de geschatte onzuiverheid. De laatste kolom toont de effectieve rekentijd in seconden. Alle simulaties werden uitgevoerd op de computercluster van de KU Leuven met p = 20 processoren. Merk opnieuw op hoe dicht de variantie van de schatter aanleunt bij de vereiste 2 /2 uit de foutenanalyse §1.1.2.

3.3

Levelafhankelijke schatters

De performantie van het MLMC-algoritme kan nog verder worden verbeterd door het invoeren van levelafhankelijke afbrekingen van de Karhunen-Loève-expansie. Deze zijn voornamelijk nuttig wanneer het stochastisch veld k(x; ω) sterk oscillerend is en varieert op een zeer kleine lengteschaal. Een ruw rooster, bijvoorbeeld ` = 0, zal dan niet in staat zijn om alle karakteristieken van k te onthullen. Door de benadering voor het stochastisch veld voortijdig af te breken kan men dit probleem vermijden. In wat volgt bestuderen we eerst het effect van het aantal termen in de KL-expansie in meer detail. Daarna bespreken we hoe men een gepaste keuze van afbrekingsstrategie moet maken. We besluiten met enkele numerieke resultaten.

3.3.1

Effect van het aantal termen in de KL-expansie

Om na te gaan wat de potentiële winst aan rekenkost is door het invoeren van levelafhankelijke schatters bestuderen we eerst het effect van het aantal termen in de KL-expansie van het stochastische veld mKL op de totale rekenkost. Figuur 3.4 toont de (genormaliseerde) rekentijd van het MLMC-algoritme met =5e-5 voor een variërend aantal termen mKL , zowel voor d = 1 (één dimensie) als d = 2 (twee dimensies). Het geval d = 1 is het modelprobleem (2.2) met λ = 0.3, σ 2 = 1, m0 = 16 en keff als eindfunctie. We vinden we een lineair verband met het aantal termen mKL . Voor het tweedimensionale geval is λ = 0.3, σ 2 = 1 en m0 = 8, en de eindfunctie is (3.5). We vinden een kwadratisch verband met het aantal termen mKL . Het 39


variantie

variantie schatter

geschatte onzuiverheid2

tijd [s]

2.521978e-01 2.478351e-02 4.999691e-07 7.029302e-08 4.441143e+01

5e-4

2.514434e-01 2.473831e-02 1.249980e-07 6.469280e-08 1.661924e+02

1e-4

2.516025e-01 2.472191e-02 4.999993e-09 3.449318e-09 4.742028e+03

5e-5

2.516914e-01 2.471553e-02 1.249998e-09 2.259572e-10 1.925337e+04

1e-3

2.530737e-01 1.548894e-02 4.997553e-07 6.354912e-09 9.048387e+00

5e-4

2.535810e-01 1.568414e-02 1.249939e-07 2.648796e-09 2.871946e+01

1e-4

2.535669e-01 1.567609e-02 4.999982e-09 3.337494e-09 6.648402e+02

5e-5

2.535926e-01 1.567762e-02 1.249995e-09 1.997243e-10 2.849176e+03

1e-3

2.479768e-01 2.296330e-02 4.999480e-07 1.759054e-07 7.436841e+01

5e-4

2.481798e-01 2.282508e-02 1.249894e-07 1.072049e-08 3.162200e+02

1e-4

2.479098e-01 2.291682e-02 4.999937e-09 7.495888e-10 8.787179e+03

5e-5

2.480226e-01 2.292926e-02 1.249998e-09 7.313016e-10 3.452724e+04

1e-3

2.487853e-01 4.630594e-03 4.997784e-07 2.234156e-09 5.205893e+00

5e-4

2.497398e-01 4.710754e-03 1.249768e-07 2.781346e-09 1.273606e+01

1e-4

2.490487e-01 4.711781e-03 4.999975e-09 1.508795e-09 2.618899e+02

5e-5

2.490292e-01 4.707611e-03 1.249989e-09 8.516310e-11 1.137104e+03

1e-3

2.575380e-01 7.844295e-02 4.999875e-07 4.046686e-07 2.377965e+02

5e-4

2.575968e-01 7.838930e-02 1.249986e-07 3.016448e-08 9.950853e+02

1e-4

2.573387e-01 7.849341e-02 4.999986e-09 1.876994e-09 2.859894e+04

5e-5

2.575596e-01 7.840566e-02 1.249999e-09 6.541216e-10 7.995231e+04

λ = 3, σ 2 = 1

λ = 0.3, σ 2 = 1

1e-3

λ = 0.3, σ 2 = 0.5 λ = 0.03, σ 2 = 1

gemiddelde waarde

λ = 0.3, σ 2 = 2

Tabel 3.3: Parameteranalyse van het modelprobleem in twee dimensies.

40

3.3. Levelafhankelijke schatters

genormaliseerde tijd

12 10

d=1 d=2

8 6 4 2 0 200

400

600

800

1000

1200

1400

mKL Figuur 3.4: Effect van het aantal termen in de KL-expansie voor d = 1 en d = 2 op de (genormaliseerde) rekentijd van het MLMC-algoritme. De tolerantie op de RMSE is =5e-5 en de parallelle implementatie met p = 20 processoren werd gebruikt. effect van het aantal termen mKL op de rekentijd van het MLMC-algoritme is dus dramatisch, en neemt toe naarmate de dimensie d = 1, 2, . . . stijgt.

3.3.2

Doelstelling

Het basisingrediënt voor multilevel Monte-Carlo is de “telescopische som” (1.6): E [QM ] = E [QM0 ] +

L X `=1

E QM` − QM`−1 = h

i

L X

E [Y` ] .

(3.6)

`=1

Zoals eerder aangehaald mag men voor Q` in (3.6) eender welke schatter gebruiken, zolang men voor het verschil Y` = Q` − Q`−1 , ` = 0 . . . L het verschil van twee dezelfde schatters voor Q` gebruikt. Men mag dus op elk level ` de coëfficiënt k(x; ω) benaderen met een verschillend aantal termen zonder daarmee de gelijkheid (3.6) te schaden. In [38] wordt aangegeven dat de performantie van de MLMC-schatter kan worden verbeterd door gebruik te maken van deze zogenaamde levelafhankelijke schatters. Het ruwe rooster (level 0) zal immers niet in staat zijn om alle kenmerken van het stochastische veld k(x, ω), dat varieert op een fijnere schaal, accuraat te benaderen. De introductie van levelafhankelijke schatters zorgt niet voor een bijkomende onzuiverheid in (3.6), maar men moet steeds oppassen om geen bijkomende modelfouten te introduceren die trager convergeren dan de discretisatiefout. Beschouwen we het Gaussiaans verdeelde stochastisch veld Z dat wordt afgebroken na mKL ∈ N termen (2.4): ZmKL (x; ω) , E[Z(x; ω)] +

m KL p X

θn fn (x)ξn (x)

(3.7)

n=1

41

3. Analyse van de MLMC-methode voor het modelprobleem Het bijhorende lognormaal verdeelde veld is dan kmKL (x; ω) , exp[ZmKL (x; ω)] (zie [38]). Beschouw dan pmKL , de oplossing van het modelprobleem − ∇ · (kmKL (x; ω)∇pmKL (x; ω)) = f (x; ω).

(3.8)

met de coëfficiënt k vervangen door de mKL -terms benadering. De convergentie van de eigenwaarden θn ten opzichte van het aantal KL-termen kan zeer traag zijn (zie Figuur 2.3). Voor een goede benadering van E[QL ] zijn dus een groot aantal termen vereist op het fijnste rooster ` = L. Wanneer de eigenwaarden {θn }n∈N geordend zijn in dalende grootte, dan zijn de bijhorende eigenfuncties geordend volgens toenemend aantal oscillaties over D. Door het veld (3.7) nu af te breken na een gegeven aantal termen wordt dus telkens de bijdrage van de meest oscillerende eigenfuncties genegeerd. Dit leidt tot zacht verlopende problemen die op ruwe roosters preciezer kunnen worden opgelost. Uit [38, Sectie 4.1] halen we dat kMω (pmKL ) − Mω (phmKL )kLp (Ω,H 1 (D)) = O(hs ), kMω (p) − Mω (pmKL )kLp (Ω,H 1 (D)) =

O(m−σ KL ),

0<s≤1

0 < σ < ∞.

Hierin is Mω ( · ) de functionaal waarin men is geïnteresseerd, pmKL de oplossing van (3.8) (de oplossing van het modelprobleem met afgebroken KL-expansie) en phmKL de gediscretiseerde oplossing van (3.8) met roosterafstand h. Om beide termen in −s

de fout kMω (p) − Mω (phmKL )k te balanceren moet mKL,` & h` σ . Voor lognormaal verdeelde stochastische velden met 1-norm covariantiefunctie is σ = s, en dus is mKL,` & h−1 ` . Men kan aantonen dat deze resultaten de asymptotische kost ( → 0) van het MLMC-algoritme niet verslechteren. Het voordeel bij het gebruik van levelafhankelijke schatters is in absolute termen en voor een vaste tolerantie . Immers, voor lognormaal verdeelde k met een onderliggende expontentiële covariantiefunctie in de 1-norm, is de optimale ruwste roosterafstand h0 ≈ λ, zie [8]. Omdat in een praktische toepassing het fijnste rooster vast ligt is dit mogelijk een beperking op het aantal levels dat men kan gebruiken. Door het invoeren van level-afhankelijke schatters is er geen beperking meer op de grootte van h0 en kan men de levelstructuur blijven gebruiken.

3.3.3

Resultaten in één dimensie

Figuur 3.5(a) toont de variantie van de eindfunctie (3.4) in één dimensie voor λ = 0.01 met een vast aantal termen in de KL-expansie. De multilevel MonteCarlo-methode is pas voordelig voor een ruwste roosterafstand kleiner dan of gelijk aan h0 = (m0 25 )−1 = 1/64. Door het toevoegen van een levelafhankelijke schatter, waarbij men als vuistregel mKL & h−1 ` hanteert, kan men het aantal levels uitbreiden, zie Figuur 3.5(b). In Figuur 3.6 werd de kostenanalyse uit Figuur 2.7 hernomen en vergeleken met een levelafhankelijke implementatie. Deze laatste levert een winst aan rekenkost met een factor 5.88 op het rooster met M` = 32 ten opzichte van standaard Monte-Carlo. De goedkoopste MLMC-schatter met vast aantal termen in de KL-expansie op dat rooster is de 2-level schatter met kost 6.1e+5. Het voordeel 42

3.4. Besluit

−5

−5

log2 V[ · ]

0

log2 V[ · ]

0

−10

−10

−15

−20

−15 Q` Q` − Q`−1 0

1

2

3

4

5

6

7

8

9

−20

Q` Q` − Q`−1 0

1

2

3

4

5

`

`

(a)

(b)

6

7

8

9

Figuur 3.5: (a) Gedrag van de variantie van Q` en Y` = Q` − Q`−1 voor λ = 0.01, σ 2 = 1, mKL = 800, m0 = 2 en N = 105 samples in één dimensie. (b) Zelfde figuur als in (a), maar met mKL = h−1 ` . van het invoeren van levelafhankelijke schatters blijkt duidelijk uit deze laatste figuur: men kan een multilevelschatter gebruiken op ruwere roosters met een lagere kost tot gevolg. Deze resultaten zijn erg belangrijk, zeker als λ → 0: uiteindelijk zal elke redelijke roosterafstand h0 als ruw worden aanzien ten opzichte van de correlatielengte λ. Door het invoeren van de levelafhankelijke afbrekingen kan de multilevelbenadering ook worden toegepast voor zeer kleine λ’s.

3.4

Besluit

Multilevel Monte-Carlo-methoden hebben het potentieel om standaard Monte-Carlosimulaties significant te versnellen. In dit hoofdstuk werden twee technieken besproken om het modelprobleem, een elliptische PDE met lognormaal verdeelde diffusiecoëfficiënt, aan te passen voor praktische doeleinden. In een eerste deel werd aangetoond dat MLMC uitstekend in parallel kan worden uitgevoerd op een multicore machine. In principe kan men perfecte schaalbaarheid verwachten: elk deterministisch sample van de elliptische PDE kan onafhankelijk van alle andere worden berekend. Dit is ook wat geobserveerd werd in enkele numerieke experimenten. In een laatste deel bespraken we de toepassing van zogenaamde level-afhankelijke schatters: door de KL-expansie voortijdig af te breken op ruwere roosters kan men de multilevelstructuur uitbreiden met ruwere roosters, ook voor zeer lage correlatielengten λ. In het middendeel werd een korte parameteranalyse van het modelprobleem uitgevoerd. Het bleek dat, zowel voor d = 1 als d = 2, de parameters α en β in Theorema 1 enkel afhangen van het probleem en de eindfunctie. Men mag dus 43


standardised cost

107

MC 2-level 3-level 4-level 5-level 2-level 3-level 4-level 5-level

106

105 101

MC MC MC MC MC MC MC MC

trunc trunc trunc trunc

102

M` Figuur 3.6: Gestandaardiseerde kost versus aantal roosterpunten M` = m` = m0 2` voor λ = 0.3, σ 2 = 1, m0 = 16 en een vaste tolerantie op de standaardafwijking van de schatter, δ =1e-3, voor schatters met vast aantal levels: (a) met mKL = 800 (volle lijn) en (b) mKL = h−1 ` (streeplijn, met label ‘trunc’). dezelfde asymptotische kost verwachten voor eender welke correlatielengte of variantie in het model voor het stochastische veld k.

44

Hoofdstuk 4

Multilevel quasi-Monte-Carlo-methoden Bij het oplossen van PDEs met stochastische coëfficiënten is men geïnteresseerd in een grootheid die kan worden voorgesteld als de verwachte waarde van een functionaal van de oplossing. In dit hoofdstuk tonen we aan dat deze grootheden zeer efficiënt kunnen worden bepaald door het combineren van (i) multilevel Monte-Carlo en (ii) quasi-Monte-Carlo (QMC). De structuur van dit hoofdstuk is als volgt: in een eerste deel situeren we het modelprobleem (2.2) in een meer algemeen kader. Daarna bespreken we rang-1-roosterpunten (rank-1 lattice points), een familie van punten voor QMC-integratie. Deze elementen worden dan gecombineerd tot een multilevel quasi-Monte-Carlo-schatter (MLQMC). Om tot een praktisch algoritme te komen voeren we een worst-case-analyse van de fout uit. Dit levert opnieuw een uitdrukking voor het optimale aantal samples op elk level. Deze uitdrukking kan dan worden gebruikt bij het formuleren van een multilevel quasi-Monte-Carlo-algoritme. Bij het algoritme hoort ook een asymptotische schatting van de bijhorende rekenkost. Het hoofdstuk wordt besloten met enkele numerieke experimenten die de superioriteit van het nieuwe algoritme illustreren.

4.1

Veralgemeende situering van het modelprobleem

In de voorgaande hoofdstukken werd de stroming van grondwater doorheen poreuze media besproken. We toonden aan dat dit modelprobleem zich herleidt tot het oplossen van een elliptische partiële differentiaalvergelijking met stochastische diffusiecoëfficiënt k(x; ω). Om de analyse van de MLQMC-schatter in §4.4.2 te veralgemenen, hernemen we hier het modelprobleem in een generieke context. Beschouw de SPDE − ∇ · (a(x, y)∇u(x, y)) = f (x)

(4.1)

in een gebied D ⊂ Rd met gegeven randvoorwaarden (e.g. u = 0 op de rand van gebied ∂D). Hierin is x ∈ D een fysische parameter (de ruimte) en y ∈ U is een aftelbaar aantal stochastische parameters (e.g. normaal verdeelde getallen in de 45

4. Multilevel quasi-Monte-Carlo-methoden KL-expansie van een stochastisch veld). De fysische dimensie van het probleem is d. Verder nemen we aan dat de diffusiecoëfficiënt a(x, y) lineair afhankelijk is van de parameters yj : a(x, y) = a ¯(x) +

∞ X

yj ψj (x).

(4.2)

j=1

De gemiddelde waarde van de diffusiecoëfficiënt noteren we als a ¯ en de functies ψj zijn de eigenfuncties uit bijvoorbeeld een KL-expansie van de covariantie-operator van het veld a(x, y). Onder milde voorwaarden op de gemiddelde waarde a ¯ en de fluctuatiecoëfficiënten yi kan men aantonen dat de som (4.2) convergeert en dat de oplossing u van (4.1) bestaat en enig is voor alle y ∈ U , zie [13]. 1 (D) → We wensen nu de verwachte waarde van een continue functionaal G : Hmix 1 R van de oplossing u(x, y) te berekenen. De ruimte Hmix is een Sobolevruimte1 met dominating mixed smoothness, zie §4.3. De gezochte verwachte waarde is dan de integraal Z I(G(u)) =

4.2

G(u(·, y))dy.

U

(4.3)

Oplossingsstrategie

De typische strategie voor het oplossen van (4.3) bestaat uit drie deelstappen: (i) Breek de oneindige som in de voorstelling van de diffusiecoëfficiënt (4.2) af na een eindig aantal termen s: a(x, y) = a ¯(x) +

s X

yj ψj (x).

j=1

Men noemt s de stochastische dimensie van het probleem. Deze waarde komt overeen met het aantal termen in de KL-expansie mKL uit Hoofdstuk 2. We noteren de verwachte waarde van dit afgebroken probleem als Is (G(us )). Dan is I(G(u)) = lim Is (G(us )). s→+∞

(ii) Discretiseer de oplossing van het afgebroken probleem us tot ush waarin h de roosterafstand in de discretisatie is. (iii) Benader de integraal Is door een kwadratuurregel Qs,N met N punten in s dimensies: Is (G(u)) ≈ Qs,N (G(ush )). (4.4) Er zijn dus drie bronnen van fouten in de gevonden oplossing: een afbrekingsfout (i), een discretisatiefout (ii) en een kwadratuurfout (iii). 1 een Sobolevruimte, vernoemd naar de Russische wiskundige Sergei Sobolev (1908-1989), is een vectorruimte van functies, uitgerust met een norm die een combinatie is van Lp -normen van de functie zelf en diens afgeleiden tot op zekere orde. De oplossingen van partiële differentiaalvergelijkingen behoren van nature tot een Sobolevruimte.

46

4.3. Roosterregels en functieruimten

4.3

Roosterregels en functieruimten

Quasi-Monte-Carlo-methoden zijn deterministische integratiemethoden die zijn ontworpen als alternatief voor Monte-Carlo-integratie. In plaats van de willekeurige punten in Monte-Carlo-simulaties neemt men een sequentie van meer uniform verdeelde punten zodat snellere convergentie kan worden bereikt. De uniformiteit wordt gemeten aan de hand van de discrepantie (4.7), daarom noemt men dergelijke puntensets ook wel lage-discrepantiepunten [14, 30]. Veronderstel dat men geïnteresseerd is in het benaderen van de integraal I gedefinieerd op een hoog-dimensionaal gebied U = [0, 1)s I=

Z U

f (y)dy.

(4.5)

Een Monte-Carlo-methode benadert deze integraal door een een kwadratuurregel met N punten en gelijke gewichten 1/N van de vorm I≈

N 1 X f (y i ). N i=1

(4.6)

De punten y i , i = 1 . . . N zijn een sequentie van s-dimensionale vectoren, onafhankelijk en uniform verdeeld in de eenheidskubus [0, 1)s . De centrale limietstelling [41] leert ons dan dat ! N cα σ 1 X f (y i ) ± √ I∈ N i=1 N met betrouwbaarheidsniveau 1−α. Het getal cα = Φ−1 (1−α/2) met Φ de cumulatieve distributiefunctie van de standaard normaalverdeling met gemiddelde 0 en variantie 1. De variantie van f (y) is σ 2 . Gegeven (een schatting voor) deze variantie σ 2 dan volgt hieruit het aantal benodigde samples om een zekere nauwkeurigheid √ te bereiken. De Monte-Carlo-schatter convergeert naar de integraal I als O(1/ N ), onafhankelijk van de dimensie s [30]. Wanneer we deze schatter (4.6) gebruiken voor de integraal (4.4) geeft dit Is (G(u)) ≈

N 1 X G(ush (·, y n )) , Qs,N (G(ush )). N n=1

Een quasi-Monte-Carlo-methode benadert (4.5) op dezelfde manier als (4.6), alleen worden de s-dimensionale integratiepunten y i niet gekozen als onafhankelijke en uniform verdeelde punten in de eenheidskubus [0, 1)s maar op een deterministische manier. Typische voorbeelden zijn Sobolpunten en digital nets [20]. Om QMC-punten te kunnen vergelijken, met elkaar en met de uniforme verdeling, voert men het begrip discrepantie in. De discrepantie van een puntenset meet de afwijking ten opzichte van de uniforme verdeling. Voor een sequentie van punten y i ∈ [0, 1)s definiëren we de lokale discrepantie discr(J, PN ) ,

#{PN ∈ J} − vol(J) N

(4.7) 47

4. Multilevel quasi-Monte-Carlo-methoden voor een subset J van U . De lokale discrepantie discr(J) is dus een maat voor het verschil tussen het werkelijke aantal punten in J en het verwachte aantal punten vol(J) in J. Voor een rechthoekige subset J(0, y 1 ) verankerd in de oorsprong definieert ∗ als [41] men de sterdiscrepantie DN ∗ DN (PN ) = sup |discr(J, PN )| J∈E ∗

met de verzameling van alle mogelijke subsets J verankerd in de oorsprong. Verschillende definities voor het begrip discrepantie zijn mogelijk, al naar gelang de gebruikte norm over de lokale discrepantie discr(J). Men noemt een puntenset PN ∗ (P ) = O(N −1 (log N )s ), met andere een lage-discrepantieset als en slechts als DN N woorden wanneer ze sneller convergeert dan standaard Monte-Carlo-punten. Figuur 4.1 vergelijkt 50 willekeurig verdeelde punten in s = 2 dimensies met evenveel punten uit een rang-1-roosterregel (4.8). De beperkende factor in de nauwkeurigheid van standaard Monte-Carlo-simulaties is het “klonteren” van de willekeurige uniform verdeelde punten uit een pseudo-random number generator of PRNG, zie [5]. De punten zijn immers onafhankelijk: ze delen geen informatie met elkaar. Het is mogelijk dat twee punten zeer dicht bij elkaar liggen en daardoor de ruimte minder uniform opvullen. Quasi-Monte-Carlo-methoden gebruiken een deterministische sequentie van punten die gecorreleerd zijn juist om die klontering te verhinderen. Merk op hoe het eenheidsvierkant in Figuur 4.1 meer uniform wordt opgevuld. In wat volgt zijn we geïnteresseerd in zogenaamde rang-1-roosterregels (rank-1 lattice rules) met punten E∗

nz y i = frac + ∆ , PN + ∆, N

(4.8)

waarin frac(·) het fractioneel deel van een getal2 is, z ∈ Zs een deterministische sdimensionale genererende vector (generating vector) en ∆ ∼ U(0, 1)s een willekeurige verschuiving (random shift) met U de uniforme verdeling. We noteren deze set van roosterpunten (lattice points) als Pn en de optelling PN + ∆ moet modulo 1 worden geïnterpreteerd. Wanneer deze set van punten wordt gebruikt voor het benaderen van Is (G(u)) (4.4) geeft dit N 1 X nz s G uh ·, frac Is (G(u)) ≈ +∆ , Qs,N (PN + ∆; G(ush )). N n=1 N

Om een foutenschatter op het resultaat te bekomen neemt men in de praktijk q verschillende sets van onafhankelijke en identiek verdeelde (iid 3 ) verschuivingen (shifts) ∆j , j = 1 . . . q: Is (G(u)) ≈ 2

1X Qs,N (PN + ∆j ; G(ush )). q j=1 q

frac(x) = x mod 1 als x > 0 en frac(x) = (1 + x) mod 1 als x < 0. Een sequentie van stochastische variabelen is onafhankelijk en identiek verdeeld (independent and identically distributed of iid) als elke stochastische variabele dezelfde kansverdeling heeft als alle andere en als alle stochastische variabelen onderling onafhankelijk zijn. 3

48

4.3. Roosterregels en functieruimten 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

1

(a)

0

0

0.2

0.4

0.6

0.8

1

(b)

Figuur 4.1: (a) 50 Monte-Carlo-punten in s = 2 dimensies (b) 50 quasi-MonteCarlo-punten in s = 2 dimensies in basis 2 met genererende vector z =[1 182667]. De verwachte waarde van de schatter is dan nog steeds de gewenste integraal Is en een foutenschatter volgt uit het betrouwbaarheidsinterval over de verschillende ∆j . Om de uniform verdeelde punten y te transformeren tot standaard normaal verdeelde punten, zoals in het modelprobleem uit de vorige hoofdstukken, kan men √ −1 gebruik maken van de inverse cumulatieve distributiefunctie: ξ = 2Φ (y i ) = i √ −1 2erf (2y i − 1) ∼ N , met N de normale verdeling. Vooraleer verder te gaan tot het construeren van een MLQMC-algoritme introduceren we enkele begrippen uit de numerieke analyse. Deze begrippen zullen van pas komen bij de foutenanalyse van de schatter. We volgen hier grotendeels [34]. De analyse van QMC integratie vereist typisch een functieruimte gedefinieerd ten opzichte van de punten y. We beschouwen de gewogen en niet-verankerde 1 1 Sobolevruimte Hmix waarin de norm van f = G(ush ) ∈ Hmix gegeven wordt door kf k2H 1 mix

=

X u⊆{1:s}

γu−1

Z 2 ∂ |u| f (y u ; y −u )dy −u dy u . |u| s−|u| ∂y [0,1] [0,1] u

Z

(4.9)

De variabelen (yj )j∈u worden voorgesteld door y u en de variabelen (yj )j∈{1:s}\u worden voorgesteld door y −u . Er is een gewichtsparameter γu geassocieerd met elke groep van variabelen y u = (yj )j∈u met indices behorende tot de verzameling u. Voor het modelprobleem zijn deze gewichten γu van het POD- of SPOD-type, zie [13] 1 and [28]. De Sobolevruimte Hmix is een speciaal geval van de Hilbertruimten die hieronder worden besproken. Definitie 4.1. (RKHS) De Hilbertruimte met reproducerende kern ( reproducing kernel Hilbert space of RKHS) H is een Hilbertruimte4 geassocieerd met een kernfunctie 4

Een Hilbertruimte, vernoemd naar de Duitse wiskundige David Hilbert (1862-1943), is de

49

4. Multilevel quasi-Monte-Carlo-methoden K die elke functie in de ruimte kan voortbrengen door middel van een inwendig product, of equivalent, elke functionaal Ty (f ) = hK(·, y, f iH in H die een puntevaluatie is in y ∈ U = [0, 1)s is begrensd (en dus continu). Met elke RKHS is een kernfunctie geassocieerd als volgt: Definitie 4.2. (voortbrengende kernfunctie) Een functie K : U × U → R is een voortbrengende kernfunctie van de Hilbertruimte H als en slechts als voor elke y ∈ U = [0, 1)s en f ∈ H geldt dat K(·, y) ∈ H

en

Ty (f ) = f (y).

Men noteert de ruimte H met voortbrengende kernfunctie K als H(K). Definitie 4.3. (verschuivingsinvariantie) Een voortbrengende kernfunctie is verschuivingsinvariant ( shift invariant) als voor elke x, y, ∆ ∈ U = [0, 1)s geldt dat K shinv (frac (x + ∆) , frac (y + ∆)) = K shinv (x, y), met K shinv de verschuivingsinvariante kernfunctie. Als een voortbrengende kernfunctie verschuivingsinvariant is, dan kan ze geschreven worden in functie van één variabele: K shinv (x, y) = K shinv (frac (x − y) , 0). Deze verschuivingsinvariante kernfunctie K shinv (x, y) is de voortbrengende kernfunctie van een verschuivingsinvariante RKHS H shinv . Nu blijkt dat het altijd mogelijk is om de RKHS H om te vormen tot een geassocieerde RKHS H shinv door de kernfunctie K uit te middelen over alle mogelijke verschuivingen: K shinv (x, y) =

Z U

K (frac (x + ∆) , frac (y + ∆)) d∆

Tot slot definiëren we de ergst-mogelijke fout (worst-case error) van de kwadratuurregel Q die de integraal I van functies f in de ruimte H(K) benadert als e2worst-case (Q, V ) ,

sup f ∈H(K) kf kH(K) ≤1

|I(f ) − Q(f )| .

(4.10)

Men kan dan aantonen dat de mean-square worst-case-error in de RKHS H(K) gebruik makende van een rang-1-roosterregel PN + ∆ met willekeurige verschuiving veralgemening van een euclidische ruimte naar een ruimte met een eindig of oneindig aantal dimensies. De ruimte is uitgerust met een norm geïnduceerd door een inwendig product. Daardoor zijn afstanden en hoeken gedefinieerd.

50

4.4. De MLQMC-schatter ∆ gelijk is aan het kwadraat van de ergst-mogelijke fout in de bijhorende ruimte H shinv met puntenset PN , zie [37]: E∆ e2worst-case (PN + ∆, K) = e2worst-case (PN , K shinv ). h

i

Bovendien geldt voor een roosterregel met q willekeurige verschuivingen dat h i 1 E∆ e2worst-case (q × PN , K) = e2worst-case (PN , K shinv ) q

waarin q × PN de q willekeurige verschuivingen van PN voorstelt, zie [37] voor details.

4.4

De MLQMC-schatter

In wat volgt bespreken we hoe de QMC-schatter kan worden gecombineerd met de multilevel Monte-Carlo-schatter tot een MLQMC-schatter.

4.4.1

Formulering van de MLQMC-schatter

We hernemen de kerngedachte van het MLMC-algoritme: beschouw een sequentie van levels M` , ` = 0, 1 . . . en noteer de oplossing van de PDE met afgebroken aantal termen in de diffusiecoëfficiënt op level ` als us` . De multilevel Monte-Carlo schatter die gebruik maakt van L levels is Is (G(u)) ≈

L X `=0

Q` G(us` ) − G(us`−1 ) , QL (G(us )).

met Q` een schatter voor het verschil G` , G(us` ) − G(us`−1 ). Verder definiëren we G(u−1 ) , 0. Als de variantie van het verschil G` kleiner is dan de variantie van G(us` ), en als samples op level ` = 0 goedkoper zijn dan op een fijner rooster, dan is de multilevelschatter goedkoper dan de standaard Monte-Carlo-schatter. Gebruiken we voor Q` een gewone Monte-Carlo-schatter, i.e. Q` (G` ) = Qs,N (G` ) dan bekomt met een multilevel Monte-Carlo-algoritme (MLMC). Als Q` benaderd wordt door een quasi-Monte-Carlo-schatter is, zoals bovenstaande roosterregel Q` (G` ) = Qs,N (PN + ∆; G` ) dan bekomt men een multilevel quasi-Monte-Carloalgoritme (MLQMC). Beschouwen we de set van verschuivingen ∆∗ = {∆0 , ∆1 . . . ∆L } dan kan men de MLQMC-schatter met gegeven rang-1-roosterregel definiëren als s QL ∗ (∆∗ ; G(u )) =

=

L X

Qs,N` (PN` + ∆` ; G` )

`=0 L X

N` 1 X nz nz G us` ·, frac − G us`−1 ·, frac + ∆` + ∆` N N N `=0 ` n=1

We beschouwen dus een verschillende verschuivingsvector op elk level `. Merk op dat us` en us`−1 opnieuw gebaseerd zijn op een sample met dezelfde willekeurige getallen. 51

4. Multilevel quasi-Monte-Carlo-methoden

4.4.2

Foutenanalyse

s De analyse van de fout I(G(u)) − QL ∗ (∆∗ ; G(u )) zal ons opnieuw een optimaal aantal samples N` op elk level ` opleveren. We schrijven

I(G(u)) −

s QL ∗ (∆∗ ; G(u ))

= I(G(u)) − = I(G(u)) − |

L X `=0 L X

Qs,N` (PN` + ∆` ; G` ) Is (G` ) +

L X `=0

`=0

{z

|

}

A

(Is − Qs,N` (PN` + ∆` ))(G` {z

}

B

= A + B(∆∗ ) waarbij we gebruik maken van de operatornotatie Q(P + ∆)(G) = Q(P + ∆; G). De eerste (deterministische) term, A, is afkomstig van de afbrekingsfout en de discretisatiefout, A = [I(G(u − usL ))] = [I(G(u − uL ))]+[I(G(uL − usL ))]. De tweede term, B, is de kwadratuurfout. De mean-square error (MSE) wordt dan 2 h i L s E∆∗ I(G(u)) − Q∗ (∆∗ ; G(u )) = E∆∗ |A + B(∆∗ )|2 h

= E∆∗ |A|2 + |B(∆∗ )|2 + 2|A||B(∆∗ )| = A2 + E∆∗ |B(∆∗ )|2 h

i

i

2

≤ [I(G(u − uL ))]2 + [I(G(uL − usL ))]

1 {z

|

+

L X `=0

|

}

E∆∗ |(Is − Qs,N` (PN` + ∆` ))(G` )|2 h

2 {z

i }

(4.11)

waarbij de verwachte waarde wordt genomen over alle verschuivingen ∆∗ . In de derde gelijkheid verdwijnt de kruisproductterm omdat de roosterregel een zuivere schatter is, i.e. E∆∗ [B(∆∗ )] = 0, de verschuivingen zijn immers iid. 1 is opnieuw afkomstig van de afbrekingsfout en de discretisatiefout, en 2 is de stochastische fout. Het is deze laatste die we wensen te minimaliseren tegen vaste kost om het optimale aantal samples N` te bepalen. In de praktijk neemt men opnieuw q verschillende sets van iid verschuivingen ∆∗,j , j = 1 . . . q: s QL q (∆∗ ; G(u )) =

q L X 1X `=0

,

L X `=0

52

q

Qs,N` (PN` + ∆`,j ; G` )

j=1

Q`∆∗ ,

(4.12)

4.4. De MLQMC-schatter waarin Q`∆∗ een schatter is voor G` gebruik makende van een gerandomiseerde rang-1-roosterregel met verschuiving ∆`,j . Er zijn dus L · q verschillende verschuivingsvectoren van lengte s. Met behulp van de ergst-mogelijke foutenanalyse uit §4.3 kan men voor deze q verschuivingen uit (4.12) de stochastische fout 2 schrijven als

2 =

L X `=0

E∆∗ |(Is − Qs,N` (PN` + ∆` ))(G` )|2 h

i

L X 1 2 1 )kG` k2H 1 ≤ eworst-case (PN` ; Hmix

≤

`=0 L X

q

mix

2 Cs,δ 1 kG` k2H 1 mix N 2δ q ` `=0

waarin een puntenset werd gebruikt met convergentie O(1/N δ ) en Cs,δ een constante die afhangt van de dimensie van de integraal s en de gebruikte puntenset PN` . Deze uitdrukking voor de stochastische fout van de MLQMC-schatter kan nu worden geminimaliseerd tegen vaste kost om het optimale aantal samples N` op elk level ` te bepalen.

4.4.3

Optimaal aantal samples

Net als in §1.1.3 kunnen we een uitdrukking vinden voor het optimale aantal samples op elk level ` zodanig dat de stochastische fout wordt geminimaliseerd. Daartoe formuleren we volgend minimalisatieprobleem: min

2 L Cs,δ 1X kG` k2H 1 mix N 2δ q `=0 `

s.t.

qN` C` = cte, ` = 0 . . . L.

N`

qN` is het totale aantal samples dat wordt genomen op level ` en C` is de kost van één sample van G` op level `. Dit is een optimalisatieprobleem met gelijkheidsbeperking dat opnieuw kan worden geschreven als een probleem zonder beperkingen. Introduceren we daartoe de Lagrangevermenigvuldiger λ dan geeft de eerste-orde optimaliteitsvoorwaarde dat −

2δ 1 C 2 kG` k2H 1 + λqC` = 0 mix q N`2δ+1 s,δ

of N` =

v u 2 kG k2 u 2δCs,δ ` H1 2δ+1 t mix

λq 2 C`

.

(4.13)

Gewone Monte-Carlo-punten convergeren als O(N −1/2 ) dus is δ = 1/2, N` = 1 en q is het aantal samples. Op die manier zijn de willekeurige punten opnieuw 53

4. Multilevel quasi-Monte-Carlo-methoden onafhankelijk. We vinden inderdaad dat het aantal samples proportioneel moet zijn p 2 kG k2 met ·/C` , zie (1.13). Daarom kiezen we Cs,δ ` H 1 h V[G` ]. We veronderstellen mix dus dat de variantie V[f ] =

Z [0,1]s

f (y) −

Z [0,1]s

f (y)dy

!2

dy

een maat is voor de norm (4.9).

4.4.4

Een praktische formule voor het aantal samples

Overeenkomstig §1.2.1 moet men (4.13) omvormen tot een praktisch bruikbaar criterium. In wat volgt nemen we aan dat QMC-punten met δ = 1 worden gebruikt. We wensen dat de kwadratuurfout 2 kleiner is dan 2 /2. Combinatie met (4.13) geeft v s u L q u1 2 V[G` ] 2 X 3 3 t 2 C Ci2 V[Gi ]. N` ≥ q 2 s,δ C `

(4.14)

i=0

Dit is ook de uitdrukking die later in de numerieke experimenten zal worden gebruikt.

4.5

Het MLQMC-algoritme

Men kan nu, analoog aan §1.2.3, een algoritme construeren gebruik makende van de resultaten (4.11), (4.12) en (4.14). Om de vereiste QMC-punten te genereren maken we gebruik van de implementatie in [36] met een genererende vector z uit [26]. Het is belangrijk om op te merken dat de roosterregel uitbreidbaar moet zijn: het totale aantal punten is a priori niet gekend, en men moet de sequentie van punten kunnen stoppen en starten bij een welbepaalde index. 1 Start met L = −1; 2 repeat 3 L ← L + 1; 4 Genereer q willekeurige verschuivingen ∆L,j ; 5 Schat V[GL ] door de steekproefvariantie voor een initieel aantal samples q × NL0 ; 6 Bereken het optimale aantal samples op elk level ` = 0 . . . L met (4.14); 7 for j = 1 . . . q do 8 Neem extra samples met verschuivingen ∆`,j op elk level ` ≤ L volgens de nieuwe N` ; 9 end 10 if L ≥ 1 then 11 Test op convergentie met (1.17); 12 end 13 until geconvergeerd; Stap (6) in het algoritme zorgt er voor dat de stochastische fout kleiner is dan 2 /2, √ en stap (11) zorgt er voor dat de bias kleiner is dan / 2. De q verschuivingen laten 54

4.6. Theoretische rekenkost van de MLQMC-schatter ons dan toe om een foutenschatter op te stellen voor het bekomen resultaat. In [12] wordt aangegeven dat in principe q, het aantal verschuivingen per level `, ook afhangt van dat level `. Er is echter maar een milde afhankelijkheid en men raadt aan om een vast aantal verschuivingen 10 < q < 30 te gebruiken. Hoe kleiner het aantal verschuivingen, hoe groter de winst aan rekenkost, maar hoe minder nauwkeurig de foutenschatter op het resultaat.

4.6

Theoretische rekenkost van de MLQMC-schatter

Net zoals in §1.1.3 kan men de asymptotische -rekenkost van de MLQMC-schatter bepalen. Deze kost is dezelfde als voor de MLMC-schatter (weliswaar met kleinere N` ) vermenigvuldigd met q, het aantal verschuivingen: s C QL q (∆∗ ; G(u )) = q

L X `=0

N` C`

Combinatie met (4.13) geeft s −1/3 C QL q (∆∗ ; G(u )) h q

L q X 3

V[GL ]C`2 .

`=0

Wanneer dus de variantie V[GL ] twee keer zo snel of sneller daalt met ` dan de kost C` stijgt ligt de dominante rekenkost op level 0. Omgekeerd, wanneer de variantie V[GL ] meer dan twee keer zo snel stijgt dan de kost C` daalt ligt de dominante kost op grootste level L. Analoog aan §1.1.3 kunnen we deze resultaten formuleren als een Theorema. Neem aan dat men over een puntenset beschikt waarvoor het verband tussen de variantie van de schatter en de variantie van de grootheid waarin men is geïnteresseerd gegeven wordt door V[Q`∆∗ ] . N`−2δ V[G` ]. Voor multilevel Monte-Carlo betekent dit dat V[Yˆ` ] . N`−1 M`−β en dus V[Y` ] . M`−β omdat V[Yˆ` ] = N`−1 V[Y` ] (zie voorwaarde (ii) in Theorema 1). Voor multilevel quasi1 Monte-Carlo met rang-1-roosterregels in Hmix veronderstellen we het verband tussen de variantie van de schatter en de variantie van de functie V[Q`∆∗ ] = N`−2 V[G` ]. Onderstaand Theorema wordt gegeven in [39], in een licht andere vorm. Voor het bewijs verwijzen we ook naar [39]. Theorema 2. Veronderstel dat er positieve constanten α, β en γ bestaan met α ≥ 12 min(β, 2δγ) en δ ∈ [1/2, ∞) en dat volgende drie voorwaarden zijn voldaan: (i) E Is (G(uhs )) − I(G(u)) . M`−α , h

i

(ii) V[Q`∆∗ ] . N`−2δ M`−β en (iii) C` . M`γ . 55

4. Multilevel quasi-Monte-Carlo-methoden Dan geldt voor elke < exp(−1) dat er een getal L en een rij {N` }L `=0 bestaan zodat de MSE e

2

s QL q (∆∗ ; G(u ))

,E

en dat C

s QL q (∆∗ ; G(u ))

.

s QL q (∆∗ ; G(u ))

 −1/δ   

− I(G(u))

−1/δ (log )1+1/(2δ)

  −1/δ−(γ−β/2)/α

2

< 2

als β > 2δγ, als β = 2δγ, als β < 2δγ.

(4.15)

Voor gewone multilevel Monte-Carlo is δ = 1/2 en dus verwachten we een kost die groeit als 1/2 als β > γ. Voor multilevel quasi-Monte-Carlo met willekeurig 1 verschoven rang-1-roosterregels in Hmix is δ = 1 en verwachten we een kost O(−1 ) als β > 2γ. In het limietgeval δ → ∞ voor zogenaamde hogere-orde QMC-methoden is het zelfs mogelijk om een exponentieel dalende kost te bereiken in Theorema 2.

4.7

Resultaten

In wat volgt bestuderen we de efficiëntie van het nieuwe MLQMC-algoritme aan aan de hand van enkele numerieke experimenten. We tonen eerst aan dat het algoritme convergeert, en dat het bovendien convergeert naar de juiste oplossing. We beschouwen het modelprobleem (2.2) in twee dimensies en modelleren de diffusiecoëfficiënt als een lognormaal verdeeld veld met correlatielengte λ = 1 en variantie σ 2 = 1 en een vast aantal termen in de KLexpansie mKL = 50. Deze waarde komt overeen met de stochastische dimensie s uit voorgaande analyse en ligt dichter bij de praktijk dan de 1400 termen uit Hoofdstuk 3. De ruwste roosterafstand is h0 = 1/8. De verdere parameters worden gespecifieerd in onderstaande tabel. parameternaam correlatielengte variantie ruwste roosterafstand aantal KL-termen aantal verschuivingen schaalfactor voor het aantal samples N` initieel aantal samples voor variantieschatting orde van convergentie van voorwaarde (i) in Theorema 2 orde van convergentie van voorwaarde (ii) in Theorema 2 toename van de rekenkost

symbool

waarde

λ σ2 h0 mKL q Cs,δ N0 α β γ

1 1 1/8 50 10 5 16 1.4956 3.8146 1.1436

De eindfunctie is de druk op x∗ = (0.5, 0.5), waarvoor de exacte waarde gekend is door onderstaande stelling, zie [21]. 56

4.7. Resultaten

10−1 12% 14% 18% 12% 14% 14% 14% 14%

10−3

10−4

10−5

10−6

10

absolute fout

absolute fout

10

−2

10−1

−2

8% 14% 16%

6%

8%

8%

2%

0%

10−3

10−4

10−5

10−3

10−2

10−6

10−3

10−2

(a) MLMC

(b) MLQMC

Figuur 4.2: Bereikte absolute fout ten opzichte van de gevraagde tolerantie op de RMS- voor n = 50 simulaties op elke tolerantie. De fout is voor elk van de gevallen normaal verdeeld (Kolmogorov-Smirnovtest). De getallen bovenaan is het percentage aan simulaties dat een fout heeft groter dan gevraagde RMS-. Stelling 4.1. Gegeven het modelprobleem (2.2) op D = [0, 1]2 dan geldt dat voor x∗ = (0.5, 0.5) E [p(x∗ ; ω)] = 0.5. Bewijs. Beschouw de bijectie γ(x) , (1 − x1 , 1 − x2 ) op D en definieer kγ (x; ω) , k(γ(x); ω) en pγ (x; ω) , 1−p(γ(x); ω). Dan voldoet het koppel (kγ , pγ ) eveneens aan het modelprobleem (2.2), immers ∇pγ (x; ω) = ∇p(γ(x); ω). Omdat het lognormaal verdeelde veld k = log Z, met Z een Gaussiaans verdeeld veld, een covariantiefunctie C(x, y) (2.3) heeft die invariant is onder de transformatie γ geldt dat kγ ∼ k en dus is E [p(x; ω)] = E [pγ (x; ω)]

= E [1 − p(γ(x); ω)]

= 1 − E [p(γ(x); ω)] . Verder geldt dat γ(x∗ ) = x∗ en dus is het gestelde bewezen. Figuur 4.2 vergelijkt de bereikte absolute fout in n = 50 realisaties van het MLMC- en MLQMC-algoritme voor verschillende toleranties op de RMSE. De fout is in alle gevallen normaal verdeeld (Kolmogorov-Smirnovtest, [31]) met gemiddelde waarde 0. Uit [32, Hoofdstuk 6] weten we dan dat 68.27% van de waarden binnen één standaardafwijking van het gemiddelde liggen. Dit blijkt inderdaad het geval te zijn, zowel voor MLMC als voor MLQMC. De getallen boven de grafiek tonen het aantal realisaties in procent met een absolute fout die toch groter is dan gevraagde RMS-. 57


107

109

5e-3 1e-3 5e-4 1e-4 5e-5 1e-5

107

105

N`

105

= = = = = =

q × N`

103

= = = = = =

5e-3 1e-3 5e-4 1e-4 5e-5 1e-5

103

101

101 0

1

2

3

4

5

6

7

0

`

1

2

3

4

5

6

7

`

(a) MLMC

(b) MLQMC

Figuur 4.3: (a) aantal samples op elk level ` voor verschillende toleranties op de RMSE voor MLMC (b) totaal aantal samples op elk level ` voor verschillende toleranties op de RMSE voor MLQMC. Merk op dat het vereiste aantal samples volgens (4.14) een ondergrens heeft gelijk aan het minimale aantal samples q × N0 voor het schatten van V[G` ]. methode

orde

MC MC* QMC QMC* MLMC MLQMC

0.38706 0.49237 0.56556 0.99856 0.49466 1.02610

Tabel 4.1: Orde van convergentie voor elke methode in Figuur 4.4. Vervolgens bestuderen we in Figuur 4.3 hoe het totale aantal samples verdeeld is over de levelstructuur. Het totale aantal samples q × N` in het MLQMC-algoritme is beduidend lager dan in het MLMC-algoritme. De ondergrens op het totale aantal samples is q × N0 , het aantal samples dat wordt genomen voor het schatten van V[Y` ]. Deze afname in aantal samples resulteert ook in een kleinere standaardkost. Figuur 4.4 vergelijkt de gestandaardiseerde rekenkost (2.17) en (2.18) voor gewone Monte-Carlo, quasi-Monte-Carlo, multilevel Monte-Carlo en multilevel quasi-MonteCarlo. De bereikte orde van convergentie staat in Tabel 4.1. Merk√op dat we voor gewone Monte-Carlo-simulatie (i.e. MLMC met L = 1) geen O(1/ cost) bereiken. Dat komt omdat de samples voor verschillende toleranties op verschillende roosters worden genomen. Immers, we wensen de RMSE te mi58

4.7. Resultaten 10−3

MC MLMC QMC MLQMC

MC* QMC*

10−4

10−5 103

104

105

106

107

108

109

1010

rekenkost Figuur 4.4: Gevraagde tolerantie ten opzichte van de gestandaardiseerde rekenkost voor elk van de vier methoden en geschatte orde van convergentie. nimaliseren, die bestaat uit een stochastische fout en een discretisatiefout. Het optimalisatiecriterium (1.11) minimaliseert deze eerste stochastische fout, onafhankelijk van de overblijvende discretisatiefout. Voor een grotere tolerantie op de RMSE kan dan een ruwer rooster worden gebruikt die de discretisatiefout toch voldoende klein houdt. Dit ruwer rooster leidt tot een kleinere standaardkost ten opzichte van dezelfde simulatie op het fijnste level L = 7. Wanneer we eenzelfde simulatie (MC*) doen waarbij de roosterafstand h constant is op level L = 7 verkrijgt men inderdaad de verwachte orde van convergentie voor een Monte-Carlo-methode. Op analoge wijze bereikt men geen O(1/cost) voor QMC, wel voor QMC* met een vast rooster. De MLMC-methode heeft een kleinere standaardkost dan QMC maar zelfde orde van convergentie. Merk op dat het resultaat voldoet aan Theorema 2: β > γ en δ = 1/2 dus is de asymptotische kost van algoritme O(−2 ). De nieuwe MLQMCmethode is een significante verbetering ten opzichte van de voorgaande methodes. Het combineren van QMC en een multilevelstructuur leidt tot een methode die wél O(1/cost) is. Er bestaat immers een α ≥ min{β/2, γ} en een β > 2γ. Uit Theorema 2 volgt dan dat de rekenkost van het MLQMC-algoritme O(−1 ) is. Het verschil in orde van convergentie tussen deze methodes komt beter tot uiting wanneer de echte rekentijden worden beschouws, zie Figuur 4.5. Door extrapolatie van de echte tijden kan men voor een gevraagde tolerantie op de RMSE de grootte-orde van de simulatieduur voor elk van de vier methoden schatten. Hoewel de MLQMC-methode convergeert naar de juiste oplossing is een foutenschatter op het resultaat niet direct toegankelijk. Voor multilevel Monte-Carlo is er een onmiddellijk verband tussen de variantie van de eindfunctie V[YL ] en de variantie van de schatter V[YˆL ], namelijk V[YˆL ] = N`−1 V[Y` ]. Voor multilevel quasi-MonteCarlo moet men de variantie van de schatter halen uit de variantie V` over de q verschuivingen van het gemiddelde van de eindfunctie G` voor N` realisaties. We vinden dat V[Q`∆∗ ] = N`−1 V` een goed criterium is voor de variantie van de schatter. 59


totale simulatietijd [ s ]

1010

10−3

95 jaar

MC QMC MLMC MLQMC

108

1 maand

106

3 dagen

104

1 uur

102

10−4

10−5

100 10−6

Figuur 4.5: Extrapolatie van de echte rekentijden voor een asymptotische tolerantie → 0. Op die manier kan men een adaptief algoritme opstellen voor Cs,δ en deze parameter verhogen totdat de variantie van de schatter kleiner is dan 2 /2. We vinden Cs,δ ≈ 5.

4.8

Besluit

In dit hoofdstuk toonden we aan dat de multilevelschatter (1.8) significant kan worden verbeterd door het combineren van MLMC en willekeurig verschoven rang-1roosterregels. Daartoe moesten we het modelprobleem (2.2) in een meer algemeen kader plaatsen en opnieuw een uitdrukking zoeken voor de fout van de resulterende schatter. De worst-case-analyse van deze fout leverde dan opnieuw een optimaal aantal samples N` . We observeerden dat de norm van de functionaal waarin men is geïnteresseerd grote overeenkomsten vertoont met de variantie. Door voor deze norm de variantie in te vullen kunnen we opnieuw een werkbaar algoritme voor MLQMC formuleren. In principe moet men immers, zoals in [29], de RMSE afschatten met de gewogen som uit een CBC-algoritme indien de functionaal behoort tot een Sobolevruimte met POD-gewichten [14, 27] of SPOD-gewichten [12, 13]. Dit leidt tot complexe uitdrukkingen voor (4.14). Bij het algoritme hoort ook een veralgemeend Theorema dat de kost voor het algoritme voorspelt, gegeven bepaalde parameters van het probleem. In numerieke experimenten bleek ten slotte duidelijk de superioriteit van het nieuwe MLQMCalgoritme ten opzichte van de bestaande methoden: MC, QMC en MLMC.

60

Hoofdstuk 5

Multilevel Monte-Carlo met meerdere eindfuncties In dit hoofdstuk wordt bekeken hoe de multilevel Monte-Carlo-methode kan worden uitgebreid voor het bepalen van meerdere eindfuncties in dezelfde simulatie. Zo’n situatie doet zich voor wanneer men bijvoorbeeld geïnteresseerd is in het drukverloop langsheen een doorsnede van het fysische integratiedomein. Het zal blijken dat de eindfunctie met de grootste multilevelcorrectieterm bepalend is voor het aantal samples op dat level. Het hoofdstuk wordt besloten met enkele numerieke experimenten waarbij het modelprobleem, de stroming van grondwater doorheen poreuze media, wordt opgelost op meer uitdagende geometrieën.

5.1

Uitbreiding naar meerdere eindfuncties

Veronderstel dat men bij het simuleren van een stochastische partiële differentiaalvergelijking niet geïnteresseerd is in de verwachte waarde van één enkele eindfunctie QM , maar in de verwachte waarde van P functionalen QM,p , p = 1 . . . P . In plaats van P afzonderlijke simulaties te doen is het mogelijk om de verwachte waarde van alle eindfuncties te bepalen in slechts één simulatie zodanig dat aan een bepaalde tolerantie p op de RMSE is voldaan. Om het optimale aantal samples op elk level te bepalen kan men opnieuw een optimalisatieprobleem oplossen. Men wenst de totale rekenkost L X `=0

N` C`

te minimaliseren waarbij de stochastische fout voor elk van de P eindfuncties aan een ongelijkheidsbeperking moet voldoen: L X `=0

N`−1 V`,p ≤

2p , p = 1 . . . P. 2

(5.1)

Hierin is V`,p = V [Q`,p − Q`−1,p ], de variantie van de multilevelcorrectieterm voor eindfunctie p, en p is de gewenste nauwkeurigheid op de RMSE voor die eindfunctie. 61

5. Multilevel Monte-Carlo met meerdere eindfuncties Dit leidt tot een optimalisatieprobleem met P beperkingen en de invoering van P Lagrangevermenigvuldigers. Het is echter niet a priori duidelijk welk van de P vermenigvuldigers actief zal zijn [18]. Een eenvoudigere en meer praktische benadering is het invoeren van V` = max p

V`,p 2p

en als enige beperking in het optimalisatieprobleem te eisen dat 1 N`−1 V` ≤ . 2 `=0

L X

Dit is opnieuw een probleem met slechts één Lagrangevermenigvuldiger en dezelfde optimale oplossing als voorheen. Concreet moet men in een simulatie dus op level ` de variantie V`,p van de eindfunctie p invullen in (1.16) waarvoor V`,p /2p maximaal is. Voor elk van de P − 1 overblijvende functionalen is (5.1) dan zeker voldaan.

5.2

Numerieke resultaten

In wat volgt worden bovenstaande resultaten toegepast op de elliptische SPDE in twee dimensies uit het modelprobleem (2.2). De meerdere eindfuncties die gezocht worden corresponderen telkens met het drukverloop p(¯ x; ω) langsheen een doorsnede ¯ van het domein. Eerst beschouwen we de elliptische PDE uit het modelprobleem x met verschillende randvoorwaarden. Daarna wordt een medium dat bestaat uit meerdere lagen bestudeerd. Tot slot bekijken we een geometrie met een centrale uitsparing.

5.2.1

Andere randvoorwaarden

Herneem de probleemstelling van de elliptische PDE −∇ · (k(x; ω)∇p(x; ω)) = f (x) op het domein [0, 1]d , d = 2 zoals in Figuur 5.1. De randvoorwaarden zijn die van een eenvoudige flow cell: een voorgeschreven druk p links en rechts en geen uitstroom op x2 = 0 en x2 = 1. k is een lognormaal verdeeld veld met correlatielengte λ = 1 en variantie σ 2 = 1. Het stochastische veld wordt voorgesteld door een KL-expansie met mKL = 1400 termen. De roosterafstand op het ruwste rooster is h0 = 1/8. Beschouw dit modelprobleem met p1 = 1 en p2 = 0. Er is geen bronterm (f ≡ 0). Figuur 5.2(a) toont de druk langsheen een doorsnede van het domein, aangeduid met de streeplijn in Figuur 5.1. De kleurcode is de (genormaliseerde) kansdichtheid van de eindfunctie zodat de maximale kans voor elke positie langsheen de doorsnede 1 is. De waarde 0 komt overeen met de 3σ-vuistregel [32]. Merk op dat de eindfunctie met maximale variantie zich in het midden van de doorsnede bevindt. Naarmate 62

5.2. Numerieke resultaten

∂p =0 ∂n 1

p1

x2

p2

0

0 x1

∂p =0 ∂n

1

Figuur 5.1: Probleemstelling van de tweedimensionale flow cell: een voorgeschreven druk links en rechts en geen uitstroom aan boven- en onderkant. men de randen nadert verkleint deze onzekerheid, om helemaal te verdwijnen op de rand. Daar is de exacte oplossing immers gegeven door de randvoorwaarden. Figuur 5.2(b) toont dezelfde resultaten voor het modelprobleem met p1 = 0 en p2 = 0. De bronterm is nu f ≡ 1. Opnieuw is de variantie van de druk p(x; ω) maximaal op x1 = 0.5 en verdwijnt de onzekerheid aan de randen. De aanwezigheid van een bronterm f verhoogt de variantie van alle eindfuncties.

5.2.2

Gelaagd medium

Een rotsformatie die vaak wordt aangetroffen in de praktijk is een gekanaliseerd medium. We verdelen D = [0, 1]2 in drie horizontale lagen en modelleren de permeabiliteit door verschillende stochastische velden. Voor de middenlaag 1/3 < x2 < 2/3 nemen we een correlatielengte λ1 = 0.1, een variantie σ12 = 1 en een gemiddelde ¯ in (2.4)) µ1 = 4. Het lognormaal verdeelde stochastisch veld waarde (i.e. log(Z) wordt gerealiseerd met behulp van een KL-expansie met 800 termen. We nemen aan dat de boven- en onderlaag uit een identiek materiaal bestaat waarvoor de permeabiliteit kan worden gemodelleerd aan de hand van een lognormaal verdeeld veld met correlatielengte λ2 = 1, variantie σ22 = 1 en gemiddelde waarde µ2 = 0. In de KL-expansie worden 512 termen gebruikt. Deze probleemstelling wordt geschetst in Figuur 5.3. Verder nemen we aan dat p1 = 0 en p2 = 0. De bronterm is f = 1. Beschouw dan een doorsnede van het domein langsheen de x2 -as op x1 = 0.5, zoals aangegeven door de streeplijn in Figuur 5.3. De druk langsheen deze doorsnede staat weergegeven in Figuur 5.4(a). De grijswaarden komen opnieuw overeen met de genormaliseerde kansdichtheidsfunctie. Merk op hoe de drielagenstructuur duidelijk zichtbaar is. In het kanaal 1/3 < x2 < 2/3 is de druk lager. De stroming kiest dan 63

5. Multilevel Monte-Carlo met meerdere eindfuncties 1.2

2

1

1.5 p(x1 )

p(x1 )

0.8 0.6 0.4

0.5

0.2 0

1

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

x1 0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

x1 0.6

0.8

(a)

1

0

0.2

0.4 (b)

Figuur 5.2: De druk langsheen een doorsnede van het domein in Figuur 5.1. (a) het geval p1 = 1, p2 = 0, f ≡ 0. (b) het geval p1 = 0, p2 = 0, f ≡ 1 deze weg met de minste weerstand, en dus kleinste druk. De variantie van de de druk in het kanaal is kleiner dan de variantie in de boven- en onderlaag. Dit resultaat is in overeenstemming met de parameteranalyse uit Hoofdstuk 3: hoe kleiner de correlatielengte, hoe kleiner de variantie van de eindfunctie. De correlatielengte λ geeft immers de afhankelijkheid weer tussen twee nabijgelegen punten van het domein. Merk op dat hier ook de druk aan de randen onzeker is, omwille van de Neumann-randvoorwaarde.

5.2.3

Andere geometrieën

Beschouw een domein zoals voorgesteld in 5.5. Dit kan een stuk rots zijn waarin in het midden een vierkant gat werd geboord met afmetingen (1/4 × 1/4). Op de randen van dit gat veronderstellen we de druk gekend, p = 1. Het lognormaal verdeelde veld heeft correlatielengte λ = 1 en variantie σ 2 = 1. De druk aan de randen p1 = 1 en p2 = 0. Er is geen bronterm. De drukverdeling langsheen een doorsnede volgens de x1 -as in x2 = 0.5 worden getoond in Figuur 5.4(b). Merk op dat de variantie van eindfunctie nul is daar waar de Dirichletvoorwaarden moeten voldaan zijn. Dit is ook het geval op de randen van de uitsparing. Om een meer diepgaande fysische analyse van het probleem te maken kan men ook de druk op verschillende punten in het ganse domein als eindfunctie kiezen. Dit laat toe om de contourlijnen van het drukverloop over het ganse gebied te visualiseren, zoals in Figuur 5.6(a). Uit dit drukprofiel kan men dan de benaderende richting van de stroomlijnen op elk van de punten bepalen. Deze zijn weergegeven door de pijlen in 64

5.2. Numerieke resultaten ∂p =0 ∂n 1

1

materiaal 1

0.8

materiaal 2

p2

6

0.6 4

x2

p1

8

0.4

materiaal 1

x2

0

0 x1

2

0.2

1

∂p =0 ∂n

0

0

0.2

0.4

0.6

0.8

0

1

x1

Figuur 5.3: Probleemstelling van de gelaagde tweedimensionale flow cell: een voorgeschreven druk links en rechts en geen uitstroom aan boven- en onderkant. De bovenste en onderste laag bestaat uit een identiek materiaal, waarvan de permeabiliteit wordt gemodelleerd door een lognormaal verdeeld stochastisch veld.

1

0.4

0.8

p(x1 )

p(x2 )

0.3

0.2

0.1

0

0.6

0.4

0.2

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

x2 0

0.2

0.4 (a)

0.6

0.8

1

0.6

0.8

1

x1 0.6

0.8

1

0

0.2

0.4 (b)

Figuur 5.4: (a) de drukverdeling langsheen een doorsnede van het domein [0, 1]2 op x1 = 0.5 volgens de x2 -richting voor het gekanaliseerde medium. (b) de drukverdeling langsheen een doorsnede van [0, 1]2 op x2 = 0.5 volgens de x1 -richting voor het probleem met centrale uitsparing.

65

5. Multilevel Monte-Carlo met meerdere eindfuncties

∂p =0 ∂n 1

p1

x2

p2

0

0 x1

∂p =0 ∂n

1

Figuur 5.5: Voorstelling van een stuk rots in het gebied [0, 1]2 met een vierkante uitsparing in het midden. De streeplijn is de doorsnede zoals getoond in Figuur 5.4(b). dezelfde figuur. De variantie van de druk op deze punten is getoond in Figuur 5.6(b). De variantie is maximaal in het punt (0.76, 0.55). Aan de randen x1 = 0, x2 = 0 en op de uitsparing is de variantie nul omwille van de Dirichletrandvoorwaarde. De rekentijd voor een dergelijke simulatie is substantieel groter (+56.2%) dan voor het geval waar slechts één eindfunctie wordt gezocht. Men moet immers op elk level en voor elk sample de oplossing interpoleren in de gevraagde punten, hier [0:0.01:1]2 . Ook de extra geheugentoegangen spelen een rol.

5.3

Besluit

In dit hoofdstuk werd aan de hand van enkele numerieke experimenten aangetoond dat de MLMC-schatter zeer goed uitbreidbaar voor het zoeken van meerdere eindfuncties met behulp van slechts één simulatie. Het is die eindfunctie met grootste variantie die bepalend is voor het aantal samples dat op elk level moet worden genomen. De enige beperkende factor is de toename aan rekentijd voor het verwerken van alle eindfuncties in elk sample op elk level. Enkele numerieke experimenten op meer uitdagende geometrieën illustreerden het gebruik van meerdere eindfuncties. Met deze uitbreiding van het MLMC-algoritme is men in staat om fysische analyses van het modelprobleem uit te voeren. Tot slot is het belangrijk om te vermelden dat men alvorens de simulatie te starten moet weten in welke eindfuncties men mogelijk geïnteresseerd is. Dit in tegenstelling tot de stochastische Galerkinmethode [10] waar men de oplossing benadert als een afgebroken polynomial chaos expansion. Één zeer omvangrijke oplossing van het bekomen stelsel levert dan een resultaat waaruit alle gewenste eindfuncties kunnen worden afgeleid. 66

5.3. Besluit

1

1

0.9

1

1

9 0.

0. 8 0. 7 0.6 0.5 0.4 0.3 0.2 0.1

0.4

1

0.4

0.6

x2

x2

1

0.6

0.2

0

0.8

0.1 0.2 0.3 0.4 0.5 0.6 0.7 .8 9 0 0.

0.8

0.2

0

0.2

0.4

0.6

x1 (a)

0.8

1

0

0

0.2

0.4

0.6

0.8

1

x1 (b)

Figuur 5.6: De contourlijnen van de verwachte waarde van de druk in discrete punten [0:0.01:1]2 (links). De pijlen geven de richting van de stroomlijnen aan. De contourlijnen van de variantie van de druk (rechts).

67

Hoofdstuk 6

Conclusies en verder werk De recente (her)ontdekking van het multilevel Monte-Carlo-algoritme is een startpunt voor veel recente onderzoeksliteratuur. In deze thesis bestudeerden we hoe het multilevelalgoritme kan worden gebruikt voor het simuleren van modellen beschreven door partiële differentiaalvergelijkingen (PDEs) met willekeurige diffusiecoëfficiënt. Het modelprobleem dat door vrijwel alle literatuur wordt gebruikt is de stroming van vloeistoffen doorheen poreuze media in een tweedimensionaal gebied. De stochastische diffusiecoëfficiënt wordt dan beschreven door een lognormaal verdeeld stochastisch veld. Omdat de kost voor het oplossen van een deterministische PDE in hogere mate afhangt van de roosterafstand dan voor stochastische differentiaalvergelijkingen is de winst aan rekenkost groter. Ook de resultaten in deze thesis bevestigen dat de multilevelstructuur zeer efficiënt is in het verlagen van de variantie van de stochastische grootheid waarin men is geïnteresseerd. Omdat een kleinere variantie zich vertaalt in minder realisaties van het stochastisch proces heeft dit een effectieve verlaging van de rekenkost tot gevolg. In Hoofdstuk 2 tonen we aan dat voor het modelprobleem in twee dimensies, met een tolerantie 1e-4 op de verwachte fout, het multilevel Monte-Carlo algoritme 75 keer sneller is dan de klassieke Monte-Carlomethode. Bovendien, hoe kleiner de gevraagde tolerantie, hoe groter de winst aan rekenkost. Gegeven een gewenst aantal realisaties van het stochastisch proces, dan kan men al deze deterministische problemen onafhankelijk van elkaar oplossen. Op die manier is het multilevelalgoritme, net als standaard Monte-Carlo, perfect parallelliseerbaar, zie Hoofdstuk 3. In dat hoofdstuk implementeren we het MLMC-algoritme op de computercluster van de KU Leuven. De parallelle implementatie van het modelprobleem in twee dimensies is tot 15 keer sneller dan de sequentiële implementatie, wanneer 20 processoren worden gebruikt. Voor het ééndimensionale geval is dit zelfs 18 keer sneller. Verder blijkt dat men uit een simulatie eveneens de kansverdeling van de gewenste stochastische grootheid kan afleiden. Alle voordelen van de standaard Monte-Carlo-schatter blijven dus behouden. In Hoofdstuk 3 tonen we aan dat de asymptotische kost voor het simuleren van een elliptische partiële differentiaalvergelijking met lognormaal verdeeld veld niet afhangt van de gebruikte parameters in dat veld, zolang de correlatielengte kleiner 69

6. Conclusies en verder werk is dan de roosterafstand op het ruwste rooster. Indien deze laatste voorwaarde niet voldaan is, kan men gebruik maken van levelafhankelijke afbrekingen in de beschrijving van het stochastische veld. In de literatuur hecht men veel belang aan de combinatie van de multilevelidee met bestaande variantiereductietechnieken, met wisselend succes. Slechts enkele veelbelovende varianten zijn Multilevel Markov chain Monte Carlo (MLMCMC), waarin observaties van de eindfunctie worden geïncorporeerd, MLMC met importance sampling [4], een variantiereductietechniek, en Continuation Multilevel Monte Carlo (CMLMC), waarbij de RMSE bestaat uit de gewogen som van de stochastische fout en de discretisatiefout [9]. In Hoofdstuk 4 combineren we de MLMC-schatter met QMC-integratie. Dit nieuwe algoritme betekent een nog grotere winst aan rekenkost, omdat de welgekozen punten zorgen voor een snellere convergentie van de verwachte waarde van de stochastische grootheid waarin men is geïnteresseerd. Voor een tolerantie 1e-4 op de verwachte fout, is het nieuwe MLQMC-algoritme nogmaals een factor 12 sneller dan de multilevelmethode. Verder illustreren we dat, voor bepaalde problemen, het mogelijk is om een rekenkost te bereiken die omgekeerd evenredig is met de gewenste tolerantie op de verwachte fout. Het formuleren van een stochastische foutenschatter ligt echter niet voor de hand. Het kan verder interessant zijn om te onderzoeken welke van deze variantiereductietechnieken moeten/kunnen worden gecombineerd om de rekenkost van het algoritme zo laag mogelijk te houden. In Hoofdstuk 5 tonen we aan dat de multilevelmethode ook toepasbaar is voor meer realistische problemen dan het modelprobleem, en dat het mogelijk is om meerdere eindfuncties met één enkele simulatie uit te rekenen. Dit wordt geïllustreerd met enkele berekeningen voor PDEs met meer algemene randvoorwaarden, met een meer complexe diffusiecoëfficiënt overeenkomstig verschillende materialen in het oplossingsdomein, en op een meer complex gebied. In elk van deze gevallen is het gedrag van de MLMC-oplossingsmethode sterk gelijkaardig aan dat van het eenvoudige modelprobleem.

70

Bĳlagen

71

Bĳlage A

Analytische uitdrukkingen voor eigenwaarden en eigenfuncties in de KL-expansie Invullen van de covariantiefunctie (2.3) voor het ééndimensionale probleem in het eigenwaardenprobleem (2.5) leidt tot de volgende integraalvergelijking Z 1 0

σ 2 e−µ|x−y| f (y)dy = θf (x),

(A.1)

met µ = 1/λ, in de literatuur vaak een Fredholm integraalvergelijking van de tweede soort genoemd. Het absolute waardeteken in (A.1) kan worden weggewerkt door het integratiegebied op te splitsen als Z x 0

2 −µ(x−y)

σ e

f (y)dy +

Z 1 x

σ 2 eµ(x−y) f (y)dy = θf (x).

(A.2)

Eénmaal afleiden naar x geeft −µ

Z x 0

2 −µ(x−y)

σ e

f (y)dy + µ

Z 1 x

σ 2 eµ(x−y) f (y)dy = θf 0 (x).

(A.3)

Merk hierbij op dat, om de afgeleide binnen het integratieteken te brengen, de integratieregel van Leibniz moet worden gebruikt: d dx

Z b(x) a(x)

g(x, t)dt =

Z b(x) a(x)

gx (x, t)dt + g(x, b(x))b0 (x) − g(x, a(x))a0 (x).

Nogmaals afleiden naar x geeft µ2

Z x 0

σ 2 e−µ(x−y) f (y)dy + µ2

Z 1 x

σ 2 eµ(x−y) f (y)dy −2µσ 2 f (x) = θf 00 (x). 73

A. Analytische uitdrukkingen voor de KL-expansie Immers, de integratieregel van Leibniz voor een tweede-orde afgeleide is d2 d2 x

Z b(x) a(x)

g(x, t)dt =

Z b(x) a(x)

gxx (x, t)dt + gx (x, b(x))b0 (x) − gx (x, a(x))a0 (x)

+ gx (x, b(x)) + gt (x, b(x))b0 (x) b0 (x) + g(x, b(x))b00 (x)

− gx (x, a(x)) + gt (x, a(x))a0 (x) a0 (x) − g(x, a(x))a00 (x)

De integraalvergelijking (A.1) kan dus worden geschreven als een tweede-orde differentiaalvergelijking µ2 θf (x) − 2µσ 2 f (x) = θf 00 (x). (A.4) Introduceer de nieuwe variabele w als

w2 =

2µσ 2 − µ2 θ θ

dan wordt de differentiaalvergelijking (A.4) f 00 (x) + w2 f (x) = 0,

0 ≤ x ≤ 1.

(A.5)

De randvoorwaarden horende bij de differentiaalvergelijking (A.5) volgen uit het evalueren van vergelijking (A.2) en (A.3) in x = 0 en x = 1. Voor x = 0 wordt dit Z 1    σ 2 e−µx2 f (y)dy 0Z 1   µ σ 2 e−µx2 f (y)dy 0

= θf (0) = θf 0 (0)

⇓

en analoog voor x = 1,

µf (0) − f 0 (0) = 0

(A.6)

µf (1) + f 0 (1) = 0

(A.7)

De integraalvergelijking (A.1) is dus omgevormd tot een gewone differentiaalvergelijking (A.5) met randvoorwaarden (A.6) en (A.7). De oplossing van de differentiaalvergelijking (A.5) is van de vorm f (x) = A sin(wx) + B cos(wx). Het toepassen van de randvoorwaarden levert volgend stelsel op: (µ tan w + w) A + (µ − w tan w) B = 0 −w A + c B = 0

(

(A.8)

Het homogene stelsel (A.8) heeft een niet-triviale oplossing als en slechts als de determinant gelijk is aan 0. Dit geeft aanleiding tot de transcendente vergelijking tan w = 74

2µw − µ2

w2

.

De oplossing van deze vergelijking is een rij wn , n = 1, 2, 3... met w1 de eerste van nul verschillende positieve oplossing. De resulterende eigenfuncties zijn dan

fn (x) = An sin(wn x) +

wn cos(wn x) µ

n = 1, 2, 3...

.

met An een normalisatieconstante zodat ||fn ||L2 [0,1] = 1. De overeenkomstige eigenwaarden zijn 2µσ 2 2λσ 2 θn = 2 = . wn + µ 2 λ2 wn2 + 1

75

Bĳlage B

De poster

77

Een parallelle Multi-Level Monte-Carlo methode voor de simulatie van stochastische partiële differentiaalvergelijkingen Probleemstelling • Uncertainty Quantification (UQ) Master Wiskundige ingenieurstechnieken

Masterproef Pieterjan Robbe

Resultaten • Situatieschets

• Stroming van grondwater doorheen poreuze media • Wet van Darcy → elliptische PDE ,

• Onzekerheid in modelparameters door gebrek aan data

• Resultaten

• Gezocht: de verwachte waarde van een functionaal van de oplossing

Promotor S. Vandewalle D. Nuyens

Figuur: domein met uitsparing

Figuur: Een typische realisatie van

Oplossing • Basis = hiërarchie van geneste roosters met roosterafstand

Academiejaar 2014-2015

• Schrijf de verwachte waarde van de eindfunctie als opeenvolgende correcties ten opzichte van het vorige ruwe rooster 10 telescoping sum”

Figuur: Verwachte waarde van de druk en stroming

Figuur: Variantie van de druk

• Vergelijking van de performantie

• Schat nu elke bijdrage in de som door gewone MC → MLMC schatter • Waarom werkt MLMC?  Samples op level  • Foutenanalyse: MSE

zijn veel goedkoper dan op level

stochastische fout aantal samples

discretisatiefout aantal levels

• MLQMC: neem integratiepunten beter dan willekeurig (rank-1 lattice rule)

Figuur: Simulatietijd voor MC, QMC, MLMC en MLQMC

Bĳlage C

Het artikel

79

1

A Practical Multilevel Quasi-Monte Carlo Method for elliptic PDEs with random coefficients Pieterjan Robbe, Dirk Nuyens and Stefan Vandewalle

Abstract—We construct a Multilevel Quasi-Monte Carlo method (MLQMC) to approximate the expected value of some functional applied to the solution of elliptic partial differential equations (PDEs) with random diffusion coefficient. The uncertainty is represented by a countably infinite number of terms. We propose an estimator for the functional of interest by combining the multilevel approach and so-called rank-1 lattice rules. The error analysis of this estimator results in an optimal number of samples at each level. We achieve a significant improvement over the classical Multilevel Monte Carlo method (MLMC). By performing a worst-case error analysis and replacing the norm of the quantity of interest by its variance we can obtain a practical implementation. Numerical results applied to the by now classical model problem of subsurface flow illustrate the superiority of the new MLQMC method over the standard MLMC method. We show that for certain problems, one can achieve a cost almost inversely proportional to the requested tolerance on the rootmean-square error.

creative in estimating the correction terms at each level. Some notable examples are level-dependent estimators [2], Markov chain Multilevel Monte Carlo [3] and Multilevel Quasi-Monte Carlo [4, 5], discussed below. The outline of the rest of this paper is as follows. Sections I-A and I-C introduce the model problem and the necessary theoretical framework. In I-B a family of QMC-points called rank-1 lattice rules is introduced. In section II we formulate the MLQMC-algorithm and present a theorem to estimate its asymptotic cost. Finally, in section III we numerically verify the gain of the MLQMC-method and show that, for certain problems, it is possible to achieve a cost almost inversely proportional to the tolerance on the rootmean-square error (RMSE).

Keywords—Uncertainty Quantification, Multilevel Quasi-Monte Carlo, stochastic partial differential equations, rank-1 lattice rules.

−∇ · (a(x, y)∇u(x, y)) = f (x) in D ⊂ Rd

I. I NTRODUCTION HE efficient numerical simulation of partial differential equations with uncertain coefficients is an important task in engineering and science. Because the uncertainty propagates through the solution, this typically leads to very high dimensional quadrature problems when computing statistics of certain quantities of interest. We focus on the problem of subsurface flow, where the uncertain diffusion coefficient is represented by a random field. Monte Carlo methods for simulating this kind of problems has been viewed as impractical due to its expense. Recently however, the Multilevel Monte Carlo method (MLMC), a combination of Monte Carlo sampling and a multigrid idea, has been successfully applied to many different problems, showing significant gains [1]. The MLMC method expresses the expected value of a certain quantity of interest on a fine grid as a telescoping sum of the expected value of that quantity on a coarse grid and a sequence of correction terms that express the difference of the quantity of interest between consecutive levels. This is possible because of the linearity of the expectation operator. Because the variance of these differences asymptotically approaches zero and samples at a coarse grid are cheaper than samples at a fine grid, one can obtain a large reduction in simulation cost. Because of the flexibility of the MLMC estimator one can be

T

Department of Computer Science, KU Leuven, Belgium, [email protected] Version of June 3, 2015.

A. Problem Formulation Consider the parametric elliptic model problem (1)

with given boundary conditions, such as u = 0 on ∂D. Here, x ∈ D is a physical parameter and y ∈ U ⊆ RN is a countably infinite number of stochastic parameters. The diffusion coefficient depends linearly on the parameters yj : a(x, y) = a ¯(x) +

∞ X

yj φj (x).

(2)

j=1

This expression arises typically from a Karhunen–Lo`eve (KL)expansion of the covariance operator of a. The mean value of the diffusion coefficient is a ¯ and φj are the eigenfunctions associated with the expansion. One can show that under mild conditions the sum (2) converges and that the solution u of (1) exists and is unique for all y ∈ U , see [6]. Our aim is now to compute the expected value of some functional 1 1 G : Hmix (D) → R of the solution u(x, y). Here, Hmix is a Sobolev space with dominating mixed smoothness, see I-C. The expected value can be expressed as the integral Z I(G(u)) = G(u(·, y))dy. (3) U

The typical solution strategy for solving (3) consists of three basic steps: (i) Truncate the infinite sum in the representation of the diffusion coefficient after s terms: s X a(x, y) = a ¯(x) + yj φj (x). j=1

We denote the expected value of this truncated problem as Is (G(us )). Hence I(G(u)) = lims→+∞ Is (G(us )).

2

(ii)

Discretize the solution of the truncated problem us to ush , where h is the grid size. (iii) Approximate the integral Is by some quadrature rule Qs,N with N points in s dimensions: Is (G(u)) ≈ Qs,N (G(ush )).

(4)

Thus, there are three sources of error: a truncation error (i), a discretization error (ii) and a quadrature error (iii). B. Rank-1 Lattice Rules Quasi-Monte Carlo (QMC) point sets are designed to serve as a deterministic alternative for the random points used in Monte Carlo methods. One chooses a sequence of more uniform distributed points to speed up convergence. Suppose we are interested in approximating the integral I defined on a high-dimensional unit cube U = [0, 1)s Z I= f (y)dy. (5) U

A quasi-Monte Carlo method approximates this integral by a quadrature rule with N points of the form I≈

N 1 X f (y i ), N i=1

(6)

see [7, 8]. Note that this is the exact same formulation as for the Monte Carlo estimator, only the points y i are no longer uniformly distributed in the unit cube [0, 1)s but are chosen as the desired deterministic sequence of QMC-points. We are interested in so-called randomly shifted rank-1 lattice rules [8], with point set nz y i = frac + ∆ , PN + ∆ N s where z ∈ Z is a deterministic s-dimensional generating vector, frac(·) is the fractional part of a number (the modulo1 operator) and ∆ ∼ U(0, 1)s is a random shift, with U the uniform distribution. We write PN for a collection of lattice points and the addition PN + ∆ is again modulo-1. Using this set of points to approximate Is (G(u)), see equation (4), yields Is (G(u)) ≈

1 N

nz G ush ·, frac +∆ N n=1 N X

, Qs,N (PN + ∆; G(ush )).

To obtain an error estimate on the result one takes q different sets of independent and identically distributed (iid) shifts ∆j , j = 1 . . . q: q

Is (G(u)) ≈

1X Qs,N (PN + ∆j ; G(ush )) q j=1

q N nz 1X 1X s G uh ·, frac + ∆j . (7) = q j=1 N n=1 N

The expected value of the estimator is still the integral Is and an error estimate follows from the confidence interval over the q sets ∆j . A generating vector z can be constructed via a component-by-component algorithm (CBC) with cost O(sN log N + s2 N ) for so-called POD weights, such that the convergence rate of the worst-case error averaged √over all shifts ∆ is close to O(1/N ), contrary to O(1/ N ) for standard Monte Carlo points. See [5] for details and a description of POD weights. Our cost model for the MLQMC estimator does not incorporate the cost of CBC. This is justified by the fact that the same lattice rule can be used for any PDE problem in a certain family. C. Function Spaces Before continuing to the formulation of a MLQMC estimator, we introduce some necessary notation and definitions used in the field of numerical analysis. These will come in handy when analysing the error of the estimator. The analysis of QMC integration typically needs a function space defined with respect to y. We consider the weighted 1 and unanchored Sobolev space Hmix with the norm of f = 1 G(ush ) ∈ Hmix given by Z 2 Z |u| X ∂ f kf k2H 1 = γu−1 (y u ; y −u )dy −u dy u mix ∂y s−|u| |u| [0,1] [0,1] u u⊆{1:s}

where y u denotes the variables (yj )j∈u and y −u are the variables (yj )j∈{1:s}\u . There is a weight parameter γu associated with each group of variables y u = (yj )j∈u with indices belonging to the set u. For our model problem these weights γu are of POD or SPOD type, see [6] and [9]. Definition 1. (RKHS) A Hilbert space H is called a reproducing kernel Hilbert space if and only if all functionals Ty (f ) = hK(·, y), f iH in H that are a point evaluation in y ∈ U = [0, 1)s are bounded, and thus continuous. The reproducing kernel associated with the Hilbert space H is a function K : U × U → R such that for every y ∈ U = [0, 1)s and f ∈ H it is true that K(·, y) ∈ H

and

Ty (f ) = f (y).

Hence, the kernel function K can reproduce any function in the Hilbert space H by the inner product with that kernel. We denote the RKHS H with reproducing kernel K as H(K). Definition 2. (Shift invariant) A reproducing kernel function is shift invariant if for every x, y, ∆ ∈ U = [0, 1)s K shinv (frac (x + ∆) , frac (y + ∆)) = K shinv (x, y),

where K shinv denotes the shift invariant kernel. If a reproducing kernel function is shift invariant, it can be written as a function of one variable: K shinv (x, y) = K shinv (frac (x − y) , 0).

This shift invariant kernel function K shinv (x, y) is the reproducing kernel function of a shift invariant RKHS H shinv . It is

3

always possible to convert the RKHS H into an associated RKHS H shinv by taking the mean over all possible shifts: Z shinv K (x, y) = K(frac (x + ∆) , frac (y + ∆))d∆. U

Finally, we define the worst-case error of some quadrature rule Q that approximates the integral I of functions f in the space H(K) as e2worst-case (Q, H) ,

sup f ∈H(K) kf kH(K) ≤1

|I(f ) − Q(f )| .

(8)

We can show that the mean-square worst-case error in the reproducing kernel Hilbert space H(K) using the rank-1 lattice rule PN + ∆ with random shift ∆ is equal to the square of the worst-case error in the corresponding shift-invariant space H shinv with point set PN : E∆ e2worst-case (PN + ∆, K) = e2worst-case (PN , K shinv ),

see [10]. Furthermore, when using q shifts, 1 E∆ e2worst-case (q × PN , K) = e2worst-case (PN , K shinv ) q where q ×PN represents the q random shifts of PN . See [10] for details. II. M ULTILEVEL Q UASI -M ONTE C ARLO METHODS Let us review the main ideas involved in the Multilevel Monte Carlo method. Consider a hierarchical sequence of grid ` sizes {h` }L `=0 , h` = h0 2 , where ` = 0 . . . L are called levels. Denote the solution of the PDE with truncated number of terms in the expansion at level ` by us` . Then, the multilevel estimator that uses L levels is L X Q` G(us` ) − G(us`−1 ) , QL (G(us )). Is (G(u)) ≈ `=0

where Q` is an unbiased estimator for the difference G` , G(us` ) − G(us`−1 ). For convenience, we define G(u−1 ) , 0. If the variance of the functional G` is smaller than the variance of G(us` ) and if samples at a coarser grid are cheaper than at a fine grid, the multilevel estimator outperforms standard Monte Carlo estimation. One can obtain a Multilevel Monte Carlo method (MLMC) by using a standard Monte Carlo estimator for Q` , i.e., Q` (G` ) = Qs,N (G` ). Similarly, a Multilevel Quasi-Monte Carlo method (MLQMC) is obtained by using a QMC estimator for Q` , such as the randomly shifted lattice rule presented above: Q` (G` ) = Qs,N (PN + ∆; G` ). Consider the set of shifts ∆∗ = {∆0 , ∆1 . . . ∆L }, a different shift vector at each level `. The MLQMC-estimator that uses L levels with shift ∆∗ is defined as s QL ∗ (∆∗ ; G(u )) =

=

L X

`=0 L X

Qs,N` (PN` + ∆` ; G` ) N` h X

nz

1 G us` ·, frac + ∆` N` n=1 N `=0 nz i + ∆` −G us`−1 ·, frac N

Note that the solutions us` and us`−1 are based on the same QMC-points. A. Error Analysis In order to find the optimal number of samples N` at each level ` that minimizes the stochastic error, we express the error s of the MLQMC estimator I(G(u)) − QL ∗ (∆∗ ; G(u )) as L X

I(G(u)) −

Qs,N` (PN` + ∆` ; G` )

`=0

= I(G(u)) −

L X

Is (G` )

`=0

L X + (Is − Qs,N` (PN` + ∆` ))(G` ) `=0

where the operator notation Q(P + ∆)(G) = Q(P + ∆; G) is used. The mean-square error (MSE) becomes h i s 2 E∆∗ I(G(u)) − QL (∆ ; G(u )) ∗ ∗ 2

2

≤ [I(G(u − uL ))] + [I(G(uL − usL ))] | {z }

1 L h i X 2 E∆∗ |(Is − Qs,N` (PN` + ∆` ))(G` )| + `=0

|

{z

2

(9)

}

where the expectation is taken over all shifts ∆∗ . The first (deterministic) term 1 is the sum of the truncation error squared and the discretization error squared. The second term,

, 2 is the stochastic error. It is the latter we want to minimize subject to a fixed cost at each level. A sufficient condition for a RMSE less than , and hence a MSE less than 2 , is that both terms in (9) are smaller than 2 /2. As for (7) we obtain an error estimate from q different sets of shifts ∆∗,j , j = 1 . . . q: s QL q (∆∗ ; G(u ))

=

q L X 1X `=0

,

L X `=0

q

Qs,N` (PN` + ∆`,j ; G` )

j=1

Q`∆∗ ,

(10)

where Q`∆∗ denotes the estimator for G` using a randomly shifted rank-1 lattice rule with shifts ∆`,j . There are thus L · q different shift vectors of length s. Using the worst-case analysis from Section I-C the stochastic error can be expressed as

2 = ≤

L X `=0

h i 2 E∆∗ |(Is − Qs,N` (PN` + ∆` ))(G` )|

L X 1 `=0

q

1 e2worst-case (PN` ; Hmix )kG` k2H 1

mix

4

−5

0

log2 V [ · ]

log2 E [| · |]

−10

−15

−5

−10

−20

−25

Q` Q` − Q`−1 0

1

2

3

−15

4

`

Q` Q` − Q`−1 0

1

2

3

4

`

Fig. 1. Estimation of the parameters β ≈ 3.7956 (left) and α ≈ 2.9352 (right). Using γ ≈ 1.5 in Theorem 1, we derive that the asymptotic cost bound is O(1/2 ) for MLMC and O(1/) for MLQMC.

≤

L X 1 `=0

q

kG` k2H 1

mix

2 Cs,δ

N`2δ

where we used a point set that has a convergence rate of O(1/N δ ) and Cs,δ is a constant which depends on the dimensionality of the integral and this point set. The computational cost of the Multilevel Quasi-Monte Carlo estimator (10) is L X s N` C` C QL (∆ ; G(u )) = q ∗ q

(11)

`=0

where C` is the cost for a single sample of G` . For the PDE application in d physical dimensions we will quantify the cost in terms of M` = h−d ` , the total number of grid points used at level `: C` = M`γ , for some γ > 0. The stochastic error of the MLQMC estimator is minimised for a fixed computational cost when s 2 kG k2 2δCs,δ ` H1 2δ+1 mix N` h , (12) 2 q C`

where “h” denotes equality up to a constant factor. For standard Multilevel Monte Carlo δ = 1/2, p N` = 1 and q is the number of samples, we find that q h ·/C` . By similarity p with the expression q h V[G` ]/C` , see e.g. [11], we propose 2 Cs,δ kGk2H 1 h V[G` ]. Combining (12) with δ = 1 and the demand that the stochastic error 2 in the multilevel estimator is less than 2 /2, we find a heuristic criterion for the number of samples that is used later in the numerical experiments: v s u 2 L q u1 2 V[G` ] X 3 2 3 t 2 N` ≥ C Ci V[Gi ]. (13) q 2 s,δ C` i=0 B. Theoretical Cost Estimate Given a point set for which V[Q`∆∗ ] = N`−2δ V[G` ] we can derive the following asymptotic cost estimates, slightly adjusted from [3, Theorem 4.1].

Theorem 1. Suppose there are positive constants α, β and γ such that α ≥ 12 min(β, 2δγ) and δ ∈ [1/2, ∞). Under the assumptions (i) E Is (G(uhs )) − I(G(u)) . M`−α , (ii) V[Q`∆∗ ] . N`−2δ M`−β and (iii) C` . M`γ there exists an L and a sequence {N` }L `=0 such that for every < exp(−1) the MSE h 2 i 2 L s s Q (∆ ; G(u )) − I(G(u)) e QL (∆ ; G(u )) , E < 2 ∗ ∗ q q and

 −1 −δ   −1 s C QL −δ (log )1+1/(2δ) q (∆∗ ; G(u )) .  −δ−1 −(γ−β/2)/α

if β > 2δγ, if β = 2δγ, if β < 2δγ. (14)

For Multilevel Monte Carlo, δ = 1/2, hence we expect a cost growing with 1/2 if β > γ. For Multilevel Quasi-Monte 1 Carlo with randomly shifted rank-1 lattice rules in Hmix ,δ=1 −1 and hence we expect a linearly growing cost O( ) if β > 2γ. If there is a point set for which δ > 1 (so-called higher-order QMC methods), it is possible to achieve exponential decaying cost in (14). C. Algorithm for MLQMC Simulation We are now ready to formulate the algorithm for MLQMCsimulation: 1) Start with L = 0 2) Generate q random shifts ∆L,j 3) Estimate V[GL ] from N 0 initial samples 4) Calculate the optimal number of samples N` , ` = 0 . . . L using (13). 5) Evaluate q × max{0, N` − N0 } samples at each level ` ≤ L using the shifts ∆`,j , j = 1 . . . q

5

TABLE I.

6) If L ≥ 1, test for convergence using Q`∆∗ h 2−α 7) If not converged, L ← L + 1 and go to step 2. Step 4 aims to minimize the stochastic error, while √ step 6 tries to ensure that the remaining bias is less than / 2. III. N UMERICAL RESULTS In this section we will apply Multilevel Quasi-Monte Carlo to elliptic PDEs with random diffusion coefficient arising in a model for subsurface flow. This steady-state single phase flow is described by Darcy’s law coupled with an incompressibility condition, leading to a (second order) PDE d

−∇ · (k(x; ω)∇p(x; ω)) = f (x; ω) in D ⊂ R .

PARAMETER VALUES USED IN THE SIMULATION

parameter name correlation length variance coarsest mesh size number of KL-modes number of shifts ∆∗,j multiplication factor for number of samples N` initial number of samples for variance estimation order of convergence of E[Q` − Q`−1 ] order of convergence of V[Q` − Q`−1 ] cost scaling factor

TABLE II.

Hence also the solution p(x; ω), the pressure gradient, will be a random field on D × Ω. A sample of k(x; ω) is produced by the (truncated) KLexpansion of Z(x; ω) , log k(x; ω), for which analytic expressions are available (see [12] or [11] for details). In practice, this infinite sum must be truncated after a certain number of modes mKL , see step (ii) in the solution strategy. The resulting deterministic PDE can be solved using any method of choice. We choose a control volume method, in which we apply a first-order finite difference method to the flux at each grid cell boundary. The resulting linear system takes the form of a five-point stencil for a 2D problem. It is sparse and ill-conditioned, but can be solved efficiently using for instance a sparse direct method. For the multilevel method we need a sequence of levels h−1 = m` , ` = 0 . . . L. Given ` a coarsest mesh size m0 we choose m` = m0 2` . The total number of grid points in a 2D problem on a square domain will then be M` = m2` . Our M ATLAB-implementation yields γ ≈ 1.5 and thus a standard cost C` = M`γ . We consider the L-shaped domain D = [0, 1]2 \[0.75, 1]2 and model the permeability k with λ = 0.3, σ = 1 and 128 terms in the KL-expansion. The source term f is sin(2πx1 ) sin(2πx2 ). Further details are specified in Table I. The quantity of interest is the effective conductivity at the left boundary x1 = 0: Z 1 ∂p Q = keff , − dx2 . (17) k ∂x1 x1 =0 0 We will first estimate the parameters α and β in Theorem 1 and next try to confirm the predicted asymptotic bound on the cost of both MLMC and MLQMC. The left plot in Figure 1 shows the behaviour of the variance of the quantity of interest Q` and the difference Q` − Q`−1 for ` = 0 . . . 4. The slope of the dashed line is approximately −3.7956, hence V[G` ] h M`−3.7956 , or β ≈ 3.7956. Similarly, the right plot shows the

value

λ σ2 h0 mKL q Cs,δ N0 α β γ

0.3 1 1/16 128 10 5 16 2.9352 3.7956 1.5

E STIMATED ORDER OF CONVERGENCE FOR FIGURE 2

(15)

Here, f (x; ω) is a source term and k(x; ω) is the diffusion coefficient. The latter is modelled as a lognormal distributed random field with covariance function kx − yk1 C(x; y) , σ 2 exp − for x, y ∈ D. (16) λ

symbol

method

cost 2.3989 1.3012 1.9632 1.1149

MC QMC MLMC MLQMC

behaviour of the expected value of Q` and Q` − Q`−1 , and we derive that α ≈ 2.9352. Since γ ≈ 1.5, we expect from Theorem 1 a cost proportional to 1/2 for MLMC and a cost proportional to 1/ for MLQMC. The only remaining unknown in the algorithm for the MLQMC simulation is the multiplication factor Cs,δ for the number of samples at each level. This constant may depend on s or δ but not on ` or N` . For the standard problem on a square domain [0, 1]2 it is known that the exact solution for a functional that is a point evaluation at x∗ = (0.5, 0.5) is 0.5, see [13]. This allows us to quantify the true error in the algorithm as a function of Cs,δ . For the experiments done with the model problem (15), we estimate Cs,δ ≈ 5. Note that an overestimation of this constant will only yield a shift in the number of samples taken, but not in the expected order of convergence of the MLQMC method. We apply standard Monte Carlo, Multilevel Monte Carlo, Quasi-Monte Carlo and Multilevel Quasi-Monte Carlo to the model problem (15). Figure 2 and Table II show the desired accuracy against the standard cost, calculated as C MLMC = N0 + for MLMC and C

MLQMC

L X

N`

`=1

= q N0 +

L X `=1

M` + M`−1 M0

M` + M`−1 N` M0

!

for MLQMC. The cost for MC is calculated as the cost for a single-level MLMC at level L and the cost for QMC as a single-level MLQMC at that same level L. First, notice how the cost for standard Monte Carlo is not quite the expected order O(2 ), since the simulation is run at different levels: a tolerance of = 0.0001 requires a 3-level estimator whereas for larger tolerances a 2-level estimator is sufficient. A smaller tolerance on the RMSE may require a

6

focus on which combination of these techniques yields the greatest cost benefit. It is for instance possible to use leveldependent truncations of the KL-expansion in estimating the quantity of interest at each level. The gain in doing so lies in absolute terms in the possibility of using even coarser mesh sizes and hence more coarser levels for a fixed tolerance.

10−2

MC QMC MLMC MLQMC

10−3

[1] [2] 10−4 103

104

105

106

107

108

109

standard cost

[3]

Fig. 2. Performance of Monte Carlo, Quasi-Monte Carlo, Multilevel Monte Carlo and Multilevel Monte Carlo: requested tolerance on RMSE versus standard cost.

[4]

finer level and hence larger standard cost. If all simulations were done at this finest grid M3 = (16 · 23 )2 = 1282 , the expected order would be exactly 2, with a larger standard cost for the larger tolerances than indicated in Figure 2. The same analysis is true for the Quasi-Monte Carlo method: we do not obtain order-1 because the simulation for smaller tolerance is run at finer levels. Note that the comparison of QMC and MLMC is not entirely fair, since QMC simulation will in practice always be performed at the finest level. Next, notice how the estimated orders for MLMC and MLQMC are in very good agreement with the theory presented above. The new MLQMC yield a significant improvement over the existing methods. Combining QMC-points with a multilevel algorithm results in a method that has a complexity nearly proportional to the requested tolerance. IV.

C ONCLUSIONS AND FURTHER WORK

In this paper we successfully applied the MLMC and MLQMC-algorithm to elliptic PDEs with random diffusion coefficient. We formulated the MLQMC-estimator using rank1 lattice rules and derived an optimal number of samples at each level ` from the error analysis. By replacing the norm of the difference Q` − Q`−1 by its variance, we obtain a practical implementation for MLQMC simulation. In principle, one must resort to theoretical bounds for the norm, as discussed in [5]. We applied this new algorithm to the by now classical problem of subsurface flow through heterogeneous materials studied in [11] and showed that the theoretical asymptotic cost bounds in Theorem 1 are valid in practice. For MLMC it is possible to achieve a cost proportional to 1/2 , where is the requested tolerance on the root-mean-square error. For MLQMC it is possible to achieve a cost bound proportional to 1/, associated with the higher order convergence rate of QMC methods. Due to the simplicity of the multilevel idea many other variance reduction techniques can be used to further reduce the cost of the MLMC estimator. Future research may

[5]

[6]

[7] [8] [9]

[10] [11]

[12] [13]

R EFERENCES M. B. Giles, “Multilevel Monte Carlo path simulation,” Operations Research, vol. 56, no. 3, pp. 607–617, 2008. A. Teckentrup, R. Scheichl, M. Giles, and E. Ullmann, “Further analysis of multilevel Monte Carlo methods for elliptic PDEs with random coefficients,” Numerische Mathematik, vol. 125, no. 3, pp. 569–600, 2013. A. L. Teckentrup, “Multilevel Monte Carlo methods and uncertainty quantification,” Ph.D. dissertation, University of Bath, 2013. [Online]. Available: http://opus.bath.ac. uk/36651/1/UnivBath PhD 2013 A Teckentrup.pdf M. B. Giles and B. J. Waterhouse, “Multilevel quasiMonte Carlo path simulation,” Advanced Financial Modelling, Radon Series on Computational and Applied Mathematics, vol. 8, pp. 165–181, 2009. F. Y. Kuo, C. Schwab, and I. H. Sloan, “Multi-level Quasi-Monte Carlo Finite Element Methods for a Class of Elliptic PDEs with Random Coefficients,” Foundations of Computational Mathematics, pp. 1–39, 2015. J. Dick, F. Y. Kuo, Q. T. Le Gia, D. Nuyens, and C. Schwab, “Higher order QMC Petrov–Galerkin Discretization for Affine Parametric Operator Equations with Random Field Inputs,” SIAM Journal on Numerical Analysis, vol. 52, no. 6, pp. 2676–2702, 2014. F. Y. Kuo and I. H. Sloan, “Lifting the curse of dimensionality,” Notices of the AMS, vol. 52, no. 11, pp. 1320– 1328, 2005. J. Dick, F. Y. Kuo, and I. H. Sloan, “High-dimensional integration: the quasi-Monte Carlo way,” Acta Numerica, vol. 22, pp. 133–288, 2013. F. Y. Kuo, C. Schwab, and I. H. Sloan, “Quasi-Monte Carlo finite element methods for a class of elliptic partial differential equations with random coefficients,” SIAM Journal on Numerical Analysis, vol. 50, no. 6, pp. 3351– 3374, 2012. D. Nuyens, “Fast construction of good lattice rules,” Ph.D. dissertation, KU Leuven, Leuven, 2007. [Online]. Available: https://lirias.kuleuven.be/handle/1979/860 K. Cliffe, M. Giles, R. Scheichl, and A. L. Teckentrup, “Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients,” Computing and Visualization in Science, vol. 14, no. 1, pp. 3–15, 2011. R. G. Ghanem and P. D. Spanos, Stochastic finite elements: a spectral approach. Courier Corporation, 2003. I. G. Graham, F. Y. Kuo, D. Nuyens, R. Scheichl, and I. H. Sloan, “Quasi-Monte Carlo methods for elliptic PDEs with random coefficients and applications,” Journal of Computational Physics, vol. 230, no. 10, pp. 3668– 3694, 2011.

Bibliografie [1] A. Barth, C. Schwab, en N. Zollinger. Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic coefficients. Numerische Mathematik, 119(1):123–161, 2011. [2] R. H. Bisseling. Parallel Scientific Computation: A structured approach using BSP and MPI. Oxford University Press, Oxford, 2004. [3] W. L. Briggs en S. F. McCormick. A multigrid tutorial. SIAM publications, Philadelphia, 2000. [4] S. Burgos en M. B. Giles. Computing Greeks using multilevel path simulation. In Monte Carlo and Quasi-Monte Carlo Methods 2010, paina’s 281–296. Springer, New York, 2012. [5] R. E. Caflisch. Monte Carlo and quasi-Monte Carlo methods. Acta Numerica, 7:1–49, 1998. [6] C. Canuto. Multilevel Monte Carlo method for PDEs with fluid dynamic applications. PhD thesis, King Abdullah University of Science and Technology, Turin, 2014. [7] J. Charrier, R. Scheichl, en A. L. Teckentrup. Finite element error analysis of elliptic PDEs with random coefficients and its application to multilevel Monte Carlo methods. SIAM Journal on Numerical Analysis, 51(1):322–352, 2013. [8] K. Cliffe, M. Giles, R. Scheichl, en A. L. Teckentrup. Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients. Computing and Visualization in Science, 14(1):3–15, 2011. [9] N. Collier, A.-L. Haji-Ali, F. Nobile, E. von Schwerin, en R. Tempone. A continuation multilevel Monte Carlo algorithm. BIT Numerical Mathematics, paina’s 1–34, 2014. [10] P. Constantine. A primer on stochastic Galerkin methods, Lecture Notes, 2007. http://inside.mines.edu/~pconstan/docs/constantine-primer.pdf. Online, geraadpleegd op 20 december 2014. 87

Bibliografie [11] B. Debusschere. Polynomial Chaos based uncertainty propagation. Slides, 4 lectures, mei 2013. KU Leuven Summer School on Uncertainty Quantification, http://people.cs.kuleuven.be/~dirk.nuyens/uqss2013/. Online, geraadpleegd op 14 september 2014. [12] J. Dick, F. Kuo, Q. T. L. Gia, en C. Schwab. Multi-level higher order QMC Galerkin discretization for affine parametric operator equations. arXiv preprint arXiv:1406.4432, 2014. [13] J. Dick, F. Y. Kuo, Q. T. Le Gia, D. Nuyens, en C. Schwab. Higher order QMC Petrov–Galerkin Discretization for Affine Parametric Operator Equations with Random Field Inputs. SIAM Journal on Numerical Analysis, 52(6):2676–2702, 2014. [14] J. Dick, F. Y. Kuo, en I. H. Sloan. High-dimensional integration: the quasi-Monte Carlo way. Acta Numerica, 22:133–288, 2013. [15] R. W. Fox, A. T. McDonald, en P. J. Pritchard. Introduction to fluid mechanics, volume 5. John Wiley & Sons, New York, 1998. [16] R. G. Ghanem en P. D. Spanos. Stochastic finite elements: a spectral approach. Courier Corporation, Chelmsford, Massachusetts, 2003. [17] M. B. Giles. Multilevel Monte Carlo Path Simulation. Operations Research, 56(3):607–617, 2008. [18] M. B. Giles. Monte Carlo Methods for Uncertainty Quantification. Slides, 4 lectures, mei 2013. KU Leuven Summer School on Uncertainty Quantification, http://people.maths.ox.ac.uk/~gilesm/mc2013/. Online, geraadpleegd op 20 augustus 2014. [19] M. B. Giles en L. Szpruch. Antithetic multilevel Monte Carlo estimation for multi-dimensional SDEs without Lévy area simulation. The Annals of Applied Probability, 24(4):1585–1620, 2014. [20] M. B. Giles en B. J. Waterhouse. Multilevel quasi-Monte Carlo path simulation. Advanced Financial Modelling, Radon Series on Computational and Applied Mathematics, 8(1):165–181, 2009. [21] I. G. Graham, F. Y. Kuo, D. Nuyens, R. Scheichl, en I. H. Sloan. Quasi-Monte Carlo methods for elliptic PDEs with random coefficients and applications. Journal of Computational Physics, 230(10):3668–3694, 2011. [22] S. Heinrich. Multilevel Monte Carlo methods. In I. Lirkov, S. D. Margenov, en J. Wasniewski, editors, Large-scale scientific computing, paina’s 58–67. Springer, Berlin, Heidelberg, 2001. [23] A. Iserles. A first course in the numerical analysis of differential equations. Cambridge University Press, Cambridge, 2009. 88

Bibliografie [24] A. Kebaier. Statistical Romberg extrapolation: a new variance reduction method and applications to option pricing. The Annals of Applied Probability, 15(4):2681–2705, 2005. [25] D. Kosambi. Statistics in function space. J. Indian Math. Soc, 7(1):76–88, 1943. [26] F. Y. Kuo. Lattice rule generating vectors, 2007. http://web.maths.unsw. edu.au/~fkuo/lattice/index.html. Online, geraadpleegd op 5 februari 2015. [27] F. Y. Kuo, C. Schwab, en I. H. Sloan. Quasi-Monte Carlo methods for highdimensional integration: the standard (weighted Hilbert space) setting and beyond. The ANZIAM Journal, 53(1):1–37, 2011. [28] F. Y. Kuo, C. Schwab, en I. H. Sloan. Quasi-Monte Carlo finite element methods for a class of elliptic partial differential equations with random coefficients. SIAM Journal on Numerical Analysis, 50(6):3351–3374, 2012. [29] F. Y. Kuo, C. Schwab, en I. H. Sloan. Multi-level Quasi-Monte Carlo Finite Element Methods for a Class of Elliptic PDEs with Random Coefficients. Foundations of Computational Mathematics, 15(2):1–39, April 2015. [30] F. Y. Kuo en I. H. Sloan. Lifting the curse of dimensionality. Notices of the AMS, 52(11):1320–1328, 2005. [31] H. W. Lilliefors. On the Kolmogorov-Smirnov test for normality with mean and variance unknown. Journal of the American Statistical Association, 62(318):399– 402, 1967. [32] P. Longley. Geographic information systems and science. John Wiley & Sons, New York, 2005. [33] N. Metropolis. The beginning of the Monte Carlo method. Los Alamos Science, 15(584):125–130, 1987. [34] R. E. Moore en M. J. Cloud. Computational Functional Analysis. Elsevier, Amsterdam, 2007. [35] K. W. Morton en D. F. Mayers. Numerical solution of partial differential equations: an introduction. Cambridge University Press, Cambridge, 2005. [36] D. Nuyens. The magic point shop of qmc point generators and generating vectors, 2010. http://people.cs.kuleuven.be/dirk.nuyens/qmc-generators/. Online, geraadpleegd op 15 december 2014. [37] D. Nuyens. Fast construction of good lattice rules. PhD thesis, KU Leuven, Leuven, 2007. https://lirias.kuleuven.be/handle/1979/860. Online, geraadpleegd op 12 februari 2014. 89

Bibliografie [38] A. Teckentrup, R. Scheichl, M. Giles, en E. Ullmann. Further analysis of multilevel Monte Carlo methods for elliptic PDEs with random coefficients. Numerische Mathematik, 125(3):569–600, 2013. [39] A. L. Teckentrup. Multilevel Monte Carlo methods and uncertainty quantification. PhD thesis, University of Bath, 2013. http://opus.bath.ac.uk/36651/1/ UnivBath_PhD_2013_A_Teckentrup.pdf. Online, geraadpleegd op 10 oktober 2014. [40] S. Torquato. Random heterogeneous materials: microstructure and macroscopic properties, volume 16. Springer Science & Business Media, New York, 2002. [41] B. Tuffin. Randomization of Quasi-Monte Carlo Methods for Error Estimation: Survey and Normal Approximation. Monte Carlo Methods and Applications, 10(3-4):617–628, 2004.

90

KU Leuven Faculteit Ingenieurswetenschappen

2014 – 2015

Fiche masterproef Student: Pieterjan Robbe Titel: Een parallelle multilevel Monte-Carlo-methode voor de simulatie van stochastische partiële differentiaalvergelijkingen Engelse titel: A parallel Multilevel Monte Carlo method for the simulation of stochastic partial differential equations UDC : 519.6 Korte inhoud: We onderzoeken hoe de multilevel Monte-Carlo-methode (MLMC) kan worden gebruik bij het simuleren van partiële differentiaalvergelijkingen met onzekerheid in de diffusiecoëfficient. Dergelijke vergelijkingen worden bijvoorbeeld gebruikt bij het modelleren van de stroming van vloeistoffen doorheen poreuze media. De onzekerheid wordt dan voorgesteld als een lognormaal verdeeld veld. We bestuderen de performantie van het algoritme en enkele varianten aan de hand van een parallelle implementatie op de computercluster van de KU Leuven. Daarnaast introduceren we een nieuw algoritme dat de multilevelschatter combineert met Quasi-Monte-Carlopunten om een snellere convergentie te bereiken. Numerieke experimenten illustreren de superioriteit van het nieuwe algoritme. We tonen aan dat het voor bepaalde problemen mogelijk is om de oplossing te berekenen, met een kost die omgekeerd evenredig is met de gevraagde tolerantie op de root-mean-square-error. Dit is heel wat sneller dan de klassieke algoritmes, waarvoor de kost omgekeerd evenredig is met de tolerantie in het kwadraat.

Thesis voorgedragen tot het behalen van de graad van Master of Science in de ingenieurswetenschappen: wiskundige ingenieurstechnieken Promotoren: Prof. Dr. ir. S. Vandewalle Prof. Dr. ir. D. Nuyens Assessoren: Prof. Dr. ir. G. Samaey Prof. Dr. R. Vandebril Begeleider: G. Suryanarayana

Een parallelle multilevel Monte-Carlo-methode voor de simulatie van stochastische partiële differentiaalvergelijkingen

Recommend Documents