[PPT] - Metropolis-Hastings Algorithms in Function Space for Bayesian PowerPoint Presentation

SLIDE 1

Metropolis-Hastings Algorithms in Function Space for Bayesian Inverse Problems

Björn Sprungk, Oliver Ernst, Hans-Jörg Starkloff, Daniel Rudolf Advances in Uncertainty Quantification Methods, Algorithms and Applications (UQAW 2015)

SRI Center for UQ in CS&E, KAUST, Saudi Arabia January 9, 2015

FAKULTÄT FÜR MATHEMATIK

SLIDE 2

Outline

1 Motivation 2 MCMC in Function Space 3 Better Proposals 4 Convergence 5 Numerical Examples

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 2 / 24

SLIDE 3

Outline

1 Motivation 2 MCMC in Function Space 3 Better Proposals 4 Convergence 5 Numerical Examples

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 3 / 24

SLIDE 4

Motivation

Uncertainty in groundwater flow modelling

6.05 6.1 6.15 6.2 x 10

5

3.57 3.575 3.58 3.585 3.59 3.595 x 10

6

UTM Easting UTM Northing

WIPP Data

transmissivity pressure head WIPP site boundary computational domain

Typical in groundwater flow modelling:

Limited observational data about

aquifer (hydraulic conductivity)

Different kinds of measurement:

direct and indirect data (transmissivity of aquifer and, e.g., pressure head)

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 4 / 24

SLIDE 5

Motivation

Uncertainty in groundwater flow modelling

0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 x 10

4

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 x 10

4

x y

Typical in groundwater flow modelling:

Limited observational data about

aquifer (hydraulic conductivity)

Different kinds of measurement:

direct and indirect data (transmissivity of aquifer and, e.g., pressure head)

Quantity of interest (QoI):

functional of flux or pressure (e.g., breakthrough time of pollutant/radionuclides)

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 4 / 24

SLIDE 6

Motivation

Uncertainty in groundwater flow modelling

4 4.5 5 5.5 6 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

log t ECDF of log travel time 20 000 MC samples

M = 30 M = 100 M = 500 M = 1000

Typical in groundwater flow modelling:

Limited observational data about

aquifer (hydraulic conductivity)

Different kinds of measurement:

direct and indirect data (transmissivity of aquifer and, e.g., pressure head)

Quantity of interest (QoI):

functional of flux or pressure (e.g., breakthrough time of pollutant/radionuclides)

UQ approach: model uncertainty

about conductivity (and therefore QoI) stochastically

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 4 / 24

SLIDE 7

Motivation

Uncertainty in groundwater flow modelling

Groundwater flow model:

−∇ · (eκ(x) ∇p(x)) = 0, p|∂D = p0

Model κ as random function from Hilbert space H , e.g., H = L2(D):

κ(x, ω) =

m≥1

ξm(ω) φm(x), {φm}∞

m=1 CONS of H .

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 5 / 24

SLIDE 8

Motivation

Uncertainty in groundwater flow modelling

Groundwater flow model:

−∇ · (eκ(x) ∇p(x)) = 0, p|∂D = p0

Model κ as random function from Hilbert space H , e.g., H = L2(D):

κ(x, ω) =

m≥1

ξm(ω) φm(x), {φm}∞

m=1 CONS of H .

Employ direct data to fit Gaussian prior model

κ ∼ µ0 = N(0, C0), resp. (ξm)m∈N = ξ ∼ µ0 = N(0, C0) on ℓ2(R).

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 5 / 24

SLIDE 9

Motivation

Uncertainty in groundwater flow modelling

Groundwater flow model:

−∇ · (eκ(x) ∇p(x)) = 0, p|∂D = p0

Model κ as random function from Hilbert space H , e.g., H = L2(D):

κ(x, ω) =

m≥1

ξm(ω) φm(x), {φm}∞

m=1 CONS of H .

Employ direct data to fit Gaussian prior model

κ ∼ µ0 = N(0, C0), resp. (ξm)m∈N = ξ ∼ µ0 = N(0, C0) on ℓ2(R).

Incorporate indirect data y by conditioning prior ξ ∼ µ0 on

y = G(κ(ξ)) + ε where G forward operator and ε ∼ N(0, Σ) Gaussian noise.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 5 / 24

SLIDE 10

Motivation

Uncertainty in groundwater flow modelling

Groundwater flow model:

−∇ · (eκ(x) ∇p(x)) = 0, p|∂D = p0

Model κ as random function from Hilbert space H , e.g., H = L2(D):

κ(x, ω) =

m≥1

ξm(ω) φm(x), {φm}∞

m=1 CONS of H .

Employ direct data to fit Gaussian prior model

κ ∼ µ0 = N(0, C0), resp. (ξm)m∈N = ξ ∼ µ0 = N(0, C0) on ℓ2(R).

Incorporate indirect data y by conditioning prior ξ ∼ µ0 on

y = G(κ(ξ)) + ε where G forward operator and ε ∼ N(0, Σ) Gaussian noise.

Bayes’ rule yields conditional/posterior distribution µ: [Stuart, 2010]

µ(dξ) ∝ e−Φ(ξ) µ0(dξ), Φ(ξ) = 1 2y − G(κ(ξ))2

Σ−1.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 5 / 24

SLIDE 11

Outline

1 Motivation 2 MCMC in Function Space 3 Better Proposals 4 Convergence 5 Numerical Examples

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 6 / 24

SLIDE 12

MCMC in Function Space

Sampling the QoI from the posterior

Markov chain (ξn)n∈N in ℓ2(R) with transition kernel Q(η, A) := P(ξn+1 ∈ A|ξn = η), A ∈ B(ℓ2(R)) which is reversible w.r.t. µ: Q(ξ, dη) µ(dξ) = Q(η, dξ) µ(dη) ⇒ µ = µQ.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 7 / 24

SLIDE 13

MCMC in Function Space

Sampling the QoI from the posterior

Markov chain (ξn)n∈N in ℓ2(R) with transition kernel Q(η, A) := P(ξn+1 ∈ A|ξn = η), A ∈ B(ℓ2(R)) which is reversible w.r.t. µ: Q(ξ, dη) µ(dξ) = Q(η, dξ) µ(dη) ⇒ µ = µQ. Then – under suitable conditions (later) – we have for QoI f 1 N

N

n=1

f (ξn)

N→∞

− − − − →

f (ξ) µ(dξ) = Eµ [f ] .

Mean squared error ∝ N− 1

2 , constant is sum of autocovariances:

∞

k=−∞

γ(k), γ(k) = Cov

f (ξ1), f (ξ1+k)
,

ξ1 ∼ µ. Rapid decay of autocovariance function γ ⇒ high statistical efficiency.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 7 / 24

SLIDE 14

MCMC in Function Space

Metropolis-Hastings MCMC

Focus on Metropolis-Hastings (MH) MCMC where ξn → ξn+1 is as follows:

1 Propose new state η according to proposal kernel q(ξn, dη), e.g.,

η ∼ q(ξn, ·) = N(ξn, s2C0), s ∈ R+ stepsize.

2 Accept proposal η with probability α(ξn, η):

draw a ∼ U[0, 1] and set ξn+1 =

η,

a ≤ α(ξn, η), ξn,

therwise.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 8 / 24

SLIDE 15

MCMC in Function Space

Metropolis-Hastings MCMC

Focus on Metropolis-Hastings (MH) MCMC where ξn → ξn+1 is as follows:

1 Propose new state η according to proposal kernel q(ξn, dη), e.g.,

η ∼ q(ξn, ·) = N(ξn, s2C0), s ∈ R+ stepsize.

2 Accept proposal η with probability α(ξn, η):

draw a ∼ U[0, 1] and set ξn+1 =

η,

a ≤ α(ξn, η), ξn,

therwise.

Then MH chain has transition kernel Q(ξ, dη) = α(ξ, η)q(ξ, dη) +

1 −
α(ξ, ζ) q(ξ, dζ)
= rejection prob.

δξ(dη),

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 8 / 24

SLIDE 16

MCMC in Function Space

MH-MCMC in Hilbert Space

Sufficient for reversibility w.r.t. µ is the choice α(ξk, η) = min

1, dν⊤

dν (ξk, η)

,

where ν(dξ, dη) := q(ξ, dη) µ(dξ), ν⊤(dξ, dη) := ν(dη, dξ).

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 9 / 24

SLIDE 17

MCMC in Function Space

MH-MCMC in Hilbert Space

Sufficient for reversibility w.r.t. µ is the choice α(ξk, η) = min

1, dν⊤

dν (ξk, η)

,

where ν(dξ, dη) := q(ξ, dη) µ(dξ), ν⊤(dξ, dη) := ν(dη, dξ). In finite dimensions dν⊤

dν is simply ratio of densities (w.r.t. Lebesgue measure).

E.g., if q has density ρ(|ξ − η|), then dν⊤ dν (ξ, η) = π(η) π(ξ) , µ(dξ) ∝ π(ξ) dξ.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 9 / 24

SLIDE 18

MCMC in Function Space

MH-MCMC in Hilbert Space

Sufficient for reversibility w.r.t. µ is the choice α(ξk, η) = min

1, dν⊤

dν (ξk, η)

,

where ν(dξ, dη) := q(ξ, dη) µ(dξ), ν⊤(dξ, dη) := ν(dη, dξ). In finite dimensions dν⊤

dν is simply ratio of densities (w.r.t. Lebesgue measure).

E.g., if q has density ρ(|ξ − η|), then dν⊤ dν (ξ, η) = π(η) π(ξ) , µ(dξ) ∝ π(ξ) dξ. In infinite dimensions, µ0-reversibility of proposal q sufficient in order that dν⊤ dν (ξ, η) exists and dν⊤ dν (ξ, η) = dµ dµ0 (η) dµ dµ0 (ξ) −1 = eΦ(ξ)−Φ(η).

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 9 / 24

SLIDE 19

MCMC in Function Space

Implications for MCMC algorithm performance

Example: 2D groundwater flow model, synthetic data. Acceptance rate vs. stepsize for increasing dimension M of ξ = (ξ1, . . . , ξM). Random walk-proposal q(ξ) = N(ξ, s2C0)

10

−3

10

−2

10

−1

10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 stepsize s Average acceptance rate M = 10 M = 40 M = 160 M = 640

average acceptance rate: ¯ α = Eν [α(ξ, η)]

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 10 / 24

SLIDE 20

MCMC in Function Space

Implications for MCMC algorithm performance

Example: 2D groundwater flow model, synthetic data. Acceptance rate vs. stepsize for increasing dimension M of ξ = (ξ1, . . . , ξM). pCN-proposal q(ξ) = N( √ 1 − s2ξ, s2C0)

10

−3

10

−2

10

−1

10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 stepsize s Average acceptance rate M = 10 M = 40 M = 160 M = 640

introduced in [Cotter et al., 2013]

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 10 / 24

SLIDE 21

Outline

1 Motivation 2 MCMC in Function Space 3 Better Proposals 4 Convergence 5 Numerical Examples

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 11 / 24

SLIDE 22

MCMC in Function Space

Implications for MCMC algorithm performance

Example: µ = N(0, C) in 2D. Different Random Walk proposals

−1.5 −1 −0.5 0.5 1 1.5 −1 −0.5 0.5 1

q(ξ) = N(ξ, s2I)

−1.5 −1 −0.5 0.5 1 1.5 −1 −0.5 0.5 1

q(ξ) = N(ξ, s2C)

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 12 / 24

SLIDE 23

MCMC in Function Space

Implications for MCMC algorithm performance

Example: µ = N(0, C) in 2D. Different Random Walk proposals

−1.5 −1 −0.5 0.5 1 1.5 −1 −0.5 0.5 1

q(ξ) = N(ξ, s2I)

−1.5 −1 −0.5 0.5 1 1.5 −1 −0.5 0.5 1

q(ξ) = N(ξ, s2C)

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 12 / 24

SLIDE 24

MCMC in Function Space

Implications for MCMC algorithm performance

Example: µ = N(0, C) in 2D. Different Random Walk proposals

50 100 150 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

lag

ACRF

X1 X2

q(ξ) = N(ξ, s2I)

50 100 150 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

ACRF

lag

X1 X2

q(ξ) = N(ξ, s2C)

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 12 / 24

SLIDE 25

MCMC in Function Space

Implications for MCMC algorithm performance

Example: µ = N(0, C) in 2D. Different Random Walk proposals

50 100 150 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

lag

ACRF

X1 X2

q(ξ) = N(ξ, s2I)

50 100 150 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

ACRF

lag

X1 X2

q(ξ) = N(ξ, s2C)

Higher statistical efficiency of proposal with same covariance as µ shown in

[Roberts & Rosenthal, 2001].

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 12 / 24

SLIDE 26

MCMC in Function Space

Gauss-Newton type approximation of posterior covariance

If forward map G : H → Rd were linear, µ0 = N(0, C0) and ε ∼ N(0, Σ), then µ = N (m, C) , C = (C −1 + G ∗Σ−1G)−1.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 13 / 24

SLIDE 27

MCMC in Function Space

Gauss-Newton type approximation of posterior covariance

If forward map G : H → Rd were linear, µ0 = N(0, C0) and ε ∼ N(0, Σ), then µ = N (m, C) , C = (C −1 + G ∗Σ−1G)−1. Idea: Gauss-Newton-type linear approximation of nonlinear G G(ξ) ≈ G(ξ) := G(ξ0) + Lξ, L = ∇G(ξ0) and use C ≈ C = (C −1 + L∗Σ−1L)−1 as proposal covariance.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 13 / 24

SLIDE 28

MCMC in Function Space

Gauss-Newton type approximation of posterior covariance

If forward map G : H → Rd were linear, µ0 = N(0, C0) and ε ∼ N(0, Σ), then µ = N (m, C) , C = (C −1 + G ∗Σ−1G)−1. Idea: Gauss-Newton-type linear approximation of nonlinear G G(ξ) ≈ G(ξ) := G(ξ0) + Lξ, L = ∇G(ξ0) and use C ≈ C = (C −1 + L∗Σ−1L)−1 as proposal covariance. Good choice for ξ0 might be the maximum a posteriori estimator: ξMAP = argminξ

Φ(ξ) + C −1/2

ξ2 .

cf. O. Ghattas’ talk.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 13 / 24

SLIDE 29

MCMC in Function Space

Posterior-informed proposals in Hilbert space

In place of prior covariance C0, use approximated posterior covariance

C = (C −1

+ Γ)−1, Γ positive, self-adjoint, bounded, otherwise arbitrary, for a Random Walk-like proposal kernel

qs(u) = N(Psu, s2

C).

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 14 / 24

SLIDE 30

MCMC in Function Space

Posterior-informed proposals in Hilbert space

In place of prior covariance C0, use approximated posterior covariance

C = (C −1

+ Γ)−1, Γ positive, self-adjoint, bounded, otherwise arbitrary, for a Random Walk-like proposal kernel

qs(u) = N(Psu, s2

C). Enforcing reversibility of kernel qs w.r.t. µ0 – as for pCN-proposal – yields Ps = C 1/2

I − s2(I + H)−1 C −1/2

, H := C 1/2 ΓC 1/2 . We call qs Gauss-Newton pCN-proposal (GNpCN).

cf. [Law, 2013], [Cui, Law & Marzouk, 2014]

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 14 / 24

SLIDE 31

MCMC in Function Space

Random Walk vs. GNpCN

Same example as for pCN, C = (C −1 + L∗Σ−1L)−1 where L = ∇G(ξMAP).

10

−3

10

−2

10

−1

10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 stepsize s Average acceptance rate M = 10 M = 40 M = 160 M = 640

Anisotropic Random Walk

qs(u) = N(u, s2

C)

10

−3

10

−2

10

−1

10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 stepsize s Average acceptance rate M = 10 M = 40 M = 160 M = 640

GNpCN

qs(u) = N(Psu, s2

C)

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 15 / 24

SLIDE 32

Outline

1 Motivation 2 MCMC in Function Space 3 Better Proposals 4 Convergence 5 Numerical Examples

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 16 / 24

SLIDE 33

Convergence

General Remarks

Markov operator Q : L2

µ(H ) → L2 µ(H ) associated with kernel Q(ξ, dη):

Qf (ξ) :=

H

f (η) Q(ξ, dη), f ∈ L2

µ(H ).

Existence of an L2-spectral gap of operator Q

0 < gap(Q) = 1 − Q − EµL2

µ→L2 µ

implies geometric ergodicity/convergence to µ in total variation norm µ − µ0QnTV ≤ C exp(−r n).

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 17 / 24

SLIDE 34

Convergence

General Remarks

Markov operator Q : L2

µ(H ) → L2 µ(H ) associated with kernel Q(ξ, dη):

Qf (ξ) :=

H

f (η) Q(ξ, dη), f ∈ L2

µ(H ).

Existence of an L2-spectral gap of operator Q

0 < gap(Q) = 1 − Q − EµL2

µ→L2 µ

implies geometric ergodicity/convergence to µ in total variation norm µ − µ0QnTV ≤ C exp(−r n).

For the pCN-proposal a (dimension-independent) L2-spectral gap was proven

under certain conditions on Φ in [Hairer et al., 2014]

Strategy: Take a comparative approach and relate spectral gap of GNpCN

to spectral gap of pCN.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 17 / 24

SLIDE 35

Convergence

Comparison Result

Theorem (Sufficient condition for spectral gap)

Let q and q be proposal kernels which admit a density ρ(ξ; η) = dq(ξ) d q(ξ)(η) and let the Markov operators Q and Q associated with the resulting MH chains be positive. If for p > 1 sup

µ(A)∈(0,1/2]

A
Ac ρp(ξ; η)

q(ξ, dη) µ(dξ) µ(A) < ∞ then gap( Q) ≥ Kp gap(Q)2p/(p−1), Kp > 0. Based on concept of conductance of Markov chain; Cheeger inequality.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 18 / 24

SLIDE 36

Convergence

Application to pCN/GNpCN

Theorem (Absolute continuity, integrability of density)

There exists a density between the pCN- and GNpCN-proposal ρs(ξ; η) := dqs(ξ) d qs(ξ)(η) which is p-integrable w.r.t. qs(ξ) for p < 3

2:

H

ρp

s (ξ; η)

qs(ξ, dη) ≤ Kp exp r(p) s2 ξ2 . Feldman-Hajek theorem, Cameron-Martin formula.

Theorem (Positivity of Markov operators)

The pCN- and GNpCN-proposal qs, qs yield positive operators Qs and Qs. Based on invariance of Gaussians w.r.t. convolution.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 19 / 24

SLIDE 37

Convergence

Main Result

Theorem (Spectral gap for restricted measure)

If the pCN-proposal yields a Markov chain with an L2-spectral gap w.r.t. µ, then so does the GNpCN-proposal w.r.t. measure µR, µR(dξ) ∝ 1BR(ξ) µ(dξ), BR := {ξ : ξ < R}, for sufficiently large R.

For sufficiently large R, µ − µRTV arbitrarily small.
Relies on (µR-reversible) pCN chain having spectral gap w.r.t. µR.
Implies convergence (to µR); rate bound for GNpCN can be obtained if

available for pCN.

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 20 / 24

SLIDE 38

Outline

1 Motivation 2 MCMC in Function Space 3 Better Proposals 4 Convergence 5 Numerical Examples

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 21 / 24

SLIDE 39

Numerical Examples

1D toy problem

1D groundwater flow model:

d dx

eκ(x) dp

dx (x)

= 0,

p(0) = 0, p(1) = 2 Data: y =

p(xj)

4

j=1 + ε,

ε ∼ N(0, σ2I), Prior: κ ∼ N(0, ∆−1)

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 22 / 24

SLIDE 40

Numerical Examples

1D toy problem

1D groundwater flow model:

d dx

eκ(x) dp

dx (x)

= 0,

p(0) = 0, p(1) = 2 Data: y =

p(xj)

4

j=1 + ε,

ε ∼ N(0, σ2I), Prior: κ ∼ N(0, ∆−1) ACRF of 3 QoIs for pCN (solid) and GNpCN (dashed)

20 40 60 80 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

σ=0.1

lag correlation QoI1 QoI2 QoI3 20 40 60 80 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

σ=0.01

lag correlation QoI1 QoI2 QoI3

QoI1 : 1

0 eκ(x) dx,

QoI2 : maxx eκ(x), QoI3 : p(0.5)

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 22 / 24

SLIDE 41

Numerical Examples

Radioactive Waste Problem (WIPP)

WIPP site: deep geological repository for radioactice waste (USA)
Model stationary groundwater flow above repository by 2D elliptic PDE.
38 measurements of transmissivity and 33 measurements of pressure head
Prior for κ: Gaussian random field with Matérn covariance function
Of interest: travel time of accidently released radionuclides

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 23 / 24

SLIDE 42

Numerical Examples

Radioactive Waste Problem (WIPP)

Importance of data and dimension

4 4.5 5 5.5 6 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 log10 travel time (years) Probability

Cumulative distribution functions: prior (solid) and posterior (dashed)

M=10 M=50 M=100 M=500 M=1000

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 23 / 24

SLIDE 43

Numerical Examples

Radioactive Waste Problem (WIPP)

Change from prior to posterior

10 10

1

10

2

10

3

10

5

10

6

10

7

10

8

10

9

m Variance

Variances of the ξm

prior post

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 23 / 24

SLIDE 44

Numerical Examples

Radioactive Waste Problem (WIPP)

ACRF in QoI for pCN (blue) and GNpCN (green)

100 200 300 400 500 600 700 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 lag correlation

ACRF for log travel time

pCN GNpCN (MAPE)

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 23 / 24

SLIDE 45

Conclusion

Dimension-independent MCMC methods needed for Bayesian inversion of

distributed quantities.

Can improve statistical efficiency of MCMC by adapting to approximate

posterior covariance.

Presented heuristic approach to capture posterior covariance and how to

make it work in Hilbert spaces.

Developed general framework to relate convergence of equivalent MCMC

methods

Convergence of GNpCN-MCMC via spectral gap, comparison theorem
To do: Extend dimension-independent algorithm to linearize at each MCMC

step.

To do: More quantitative comparison of gap sizes

Oliver Ernst (TU Chemnitz) Metropolis-Hastings Algorithms KAUST, January 2015 24 / 24