Transport methods for sampling: low-dimensional structure and - - PowerPoint PPT Presentation

Transport methods for sampling: low-dimensional structure and preconditioning

Youssef Marzouk joint work with Daniele Bigoni, Matthew Parno, Alessio Spantini, & Olivier Zahm

Department of Aeronautics and Astronautics Center for Computational Engineering Statistics and Data Science Center Massachusetts Institute of Technology http://uqgroup.mit.edu

Support from AFOSR, DARPA, DOE

12 July 2019

Marzouk et al. MIT 1 / 40

Motivation: Bayesian inference in large-scale models

Observations y Parameters x

πpos(x) := π(x|y) ∝ π(y|x)πpr(x)

• Bayes’ rule

◮ Characterize the posterior distribution (density πpos) ◮ This is a challenging task since:

◮ x ∈ Rn is typically high-dimensional (e.g., a discretized function) ◮ πpos is non-Gaussian ◮ evaluations of the likelihood (hence πpos) may be expensive

◮ πpos can be evaluated up to a normalizing constant

Marzouk et al. MIT 2 / 40

Computational challenges

◮ Extract information from the posterior (means, covariances, event

probabilities, predictions) by evaluating posterior expectations: Eπpos[h(x)] =

• h(x)πpos(x)dx

◮ Key strategy for making this computationally tractable:

◮ Efficient and structure-exploiting sampling schemes Marzouk et al. MIT 3 / 40

Computational challenges

◮ Extract information from the posterior (means, covariances, event

probabilities, predictions) by evaluating posterior expectations: Eπpos[h(x)] =

• h(x)πpos(x)dx

◮ Key strategy for making this computationally tractable:

◮ Efficient and structure-exploiting sampling schemes ◮ This talk: relate to notions of coupling and transport. . . Marzouk et al. MIT 3 / 40

Deterministic couplings of probability measures

T η π

Core idea

◮ Choose a reference distribution η (e.g., standard Gaussian) ◮ Seek a transport map T : Rn → Rn such that T♯η = π

Marzouk et al. MIT 4 / 40

Deterministic couplings of probability measures

S =T −1 η π

Core idea

◮ Choose a reference distribution η (e.g., standard Gaussian) ◮ Seek a transport map T : Rn → Rn such that T♯η = π ◮ Equivalently, find S = T −1 such that S♯π = η

Marzouk et al. MIT 4 / 40

Deterministic couplings of probability measures

T η π

Core idea

◮ Choose a reference distribution η (e.g., standard Gaussian) ◮ Seek a transport map T : Rn → Rn such that T♯η = π ◮ Equivalently, find S = T −1 such that S♯π = η ◮ In principle, enables exact (independent, unweighted) sampling!

Marzouk et al. MIT 4 / 40

Deterministic couplings of probability measures

T η π

Core idea

◮ Choose a reference distribution η (e.g., standard Gaussian) ◮ Seek a transport map T : Rn → Rn such that T♯η = π ◮ Equivalently, find S = T −1 such that S♯π = η ◮ Satisfying these conditions only approximately can still be useful!

Marzouk et al. MIT 4 / 40

Choice of transport map

A useful building block is the Knothe–Rosenblatt rearrangement: T(x) =      T 1(x1) T 2(x1, x2) . . . T n(x1, x2, . . . , xn)     

◮ Unique triangular and monotone map satisfying T♯η = π for

absolutely continuous η, π on Rn

◮ Jacobian determinant easy to evaluate ◮ Monotonicity is essentially one-dimensional: ∂xkT k > 0 ◮ “Exposes” marginals, enables conditional sampling. . .

Marzouk et al. MIT 5 / 40

Choice of transport map

A useful building block is the Knothe–Rosenblatt rearrangement: T(x) =      T 1(x1) T 2(x1, x2) . . . T n(x1, x2, . . . , xn)     

◮ Unique triangular and monotone map satisfying T♯η = π for

absolutely continuous η, π on Rn

◮ Jacobian determinant easy to evaluate ◮ Monotonicity is essentially one-dimensional: ∂xkT k > 0 ◮ “Exposes” marginals, enables conditional sampling. . . ◮ Numerical approximations can employ a monotone parameterization

guaranteeing ∂xkT k > 0. For example: T k(x1, . . . , xk) = ak(x1, . . . , xk−1)+ xk exp(bk(x1, . . . , xk−1, w)) dw

Marzouk et al. MIT 5 / 40

How to construct triangular maps?

Construction #1: “maps from densities,” i.e., variational characterization of the direct map T [Moselhy & M 2012]

Marzouk et al. MIT 6 / 40

How to construct triangular maps?

Construction #1: “maps from densities,” i.e., variational characterization of the direct map T [Moselhy & M 2012] min

T∈T h

△

DKL( T♯ η || π ) = min

T∈T h

△

DKL( η || T −1

♯

π )

◮ π is the “target” density on Rn; η is, e.g., N(0, In) ◮ T h △ is a set of monotone lower triangular maps

◮ T h→∞

△

contains the Knothe–Rosenblatt rearrangement

◮ Expectation is with respect to the reference measure η

◮ Compute via, e.g., Monte Carlo, sparse quadrature

◮ Use unnormalized evaluations of π and its gradients ◮ No MCMC or importance sampling ◮ In general non-convex, unless π is log-concave

Marzouk et al. MIT 6 / 40

Illustrative example

min

T Eη[ − log π ◦ T −

• k

log ∂xkT k ]

◮ Parameterized map T ∈ T h △ ⊂ T△ ◮ Optimize over coefficients of

parameterization

◮ Use gradient-based optimization ◮ The posterior is in the tail of the reference

3 3 3 3 6 9 12 15 18

Marzouk et al. MIT 6 / 40

Illustrative example

min

T Eη[ − log π ◦ T −

• k

log ∂xkT k ]

◮ Parameterized map T ∈ T h △ ⊂ T△ ◮ Optimize over coefficients of

parameterization

◮ Use gradient-based optimization ◮ The posterior is in the tail of the reference

3 3 3 3 6 9 12 15 18

Marzouk et al. MIT 6 / 40

Illustrative example

min

T Eη[ − log π ◦ T −

• k

log ∂xkT k ]

◮ Parameterized map T ∈ T h △ ⊂ T△ ◮ Optimize over coefficients of

parameterization

◮ Use gradient-based optimization ◮ The posterior is in the tail of the reference

3 3 3 3 6 9 12 15 18

Marzouk et al. MIT 6 / 40

Illustrative example

min

T Eη[ − log π ◦ T −

• k

log ∂xkT k ]

◮ Parameterized map T ∈ T h △ ⊂ T△ ◮ Optimize over coefficients of

parameterization

◮ Use gradient-based optimization ◮ The posterior is in the tail of the reference

3 3 3 3 6 9 12 15 18

Marzouk et al. MIT 6 / 40

Useful features

◮ Move samples; don’t just reweigh them ◮ Independent and cheap samples: xi ∼ η ⇒ T(xi) ◮ Clear convergence criterion, even with unnormalized target density:

DKL( T♯ η || π ) ≈ 1 2 Varη

• log

η T −1

♯

¯ π

• Marzouk et al.

MIT 7 / 40

Useful features

◮ Move samples; don’t just reweigh them ◮ Independent and cheap samples: xi ∼ η ⇒ T(xi) ◮ Clear convergence criterion, even with unnormalized target density:

DKL( T♯ η || π ) ≈ 1 2 Varη

• log

η T −1

♯

¯ π

• ◮ Can either accept bias or reduce it by:

◮ Increasing the complexity of the map T ∈ T h

△

◮ Sampling the pullback T −1

♯

π using MCMC or importance sampling (more on this later)

Marzouk et al. MIT 7 / 40

Useful features

◮ Move samples; don’t just reweigh them ◮ Independent and cheap samples: xi ∼ η ⇒ T(xi) ◮ Clear convergence criterion, even with unnormalized target density:

DKL( T♯ η || π ) ≈ 1 2 Varη

• log

η T −1

♯

¯ π

• ◮ Can either accept bias or reduce it by:

◮ Increasing the complexity of the map T ∈ T h

△

◮ Sampling the pullback T −1

♯

π using MCMC or importance sampling (more on this later)

◮ Related transport constructions for inference and sampling: Stein variational

gradient descent [Liu & Wang 2016, DeTommaso 2018], normalizing flows [Rezende & Mohamed 2015], SOS polynomial flow [Jaini et al. 2019], Gibbs flow [Heng et al. 2015], particle flow filter [Reich 2011], implicit sampling [Chorin et al. 2009–2015], etc.

Marzouk et al. MIT 7 / 40

How to construct triangular maps?

Construction #2: “maps from samples” min

S∈Sh

△

DKL( S♯π || η ) = min

S∈Sh

△

DKL( π || S−1

♯

η )

◮ Suppose we have Monte Carlo samples {xi}M i=1 ∼ π ◮ For standard Gaussian η, this problem is convex and separable ◮ This is density estimation via transport! (cf. Tabak & Turner 2013)

Marzouk et al. MIT 8 / 40

How to construct triangular maps?

Construction #2: “maps from samples” min

S∈Sh

△

DKL( S♯π || η ) = min

S∈Sh

△

DKL( π || S−1

♯

η )

◮ Suppose we have Monte Carlo samples {xi}M i=1 ∼ π ◮ For standard Gaussian η, this problem is convex and separable ◮ This is density estimation via transport! (cf. Tabak & Turner 2013) ◮ Equivalent to maximum likelihood estimation of S

• S ∈ arg max

S∈Sh

△

1 M

M

• i=1

log S−1

♯

η

pullback

(xi), η = N(0, In),

◮ Each component

Sk of S can be computed separately, via smooth convex optimization

• Sk ∈ arg

min

Sk∈Sh

△,k

1 M

M

• i=1

1 2Sk(xi)2 − log ∂kSk(xi)

• Marzouk et al.

MIT 8 / 40

Low-dimensional structure of transport maps

Underlying challenge: maps in high dimensions

◮ Major bottleneck: representation of the map, e.g., cardinality of the

map basis

◮ How to make the construction/representation of high-dimensional

transports tractable?

Marzouk et al. MIT 9 / 40

Low-dimensional structure of transport maps

Underlying challenge: maps in high dimensions

◮ Major bottleneck: representation of the map, e.g., cardinality of the

map basis

◮ How to make the construction/representation of high-dimensional

transports tractable? Main ideas:

1

Exploit Markov structure of the target distribution

◮ Leads to sparsity and/or decomposability of transport maps [Spantini,

Bigoni, & M JMLR 2018]

2

Exploit certain low rank structure

◮ Near-identity or “lazy” maps [Bigoni et al. arXiv:1906.00031] Marzouk et al. MIT 9 / 40

Markov random fields

◮ Let Z1, . . . , Zn be random variables with joint density π > 0

A B S

(i, j) / ∈ E iff Zi ⊥ ⊥ Zj | ZV\{i,j}

◮ G = (V, E) encodes conditional independence (an I-map for π) ◮ Theorem [SBM 2018]: Define G s.t. (i, j) /

∈ E iff ∂xi,xj log π = 0. Then the resulting G is the unique minimal I-map for π.

Marzouk et al. MIT 10 / 40

Sparsity of transport maps

◮ Focus on the inverse triangular map S, where S♯π = η ◮ Theorem [SBM 2018]: S (a nonlinear function) inherits the same

sparsity pattern as the Cholesky factor of the incidence matrix (properly scaled) of a graphical model for π, provided that η(x) =

i η(xi)

S(x) =        S1(x1) S2(x1, x2) S3(x1, x2, x3) . . . Sn(x1, x2, . . . , xn)        = ⇒        S1(x1) S2(x1, x2) S3(x1,x2, x3) . . . Sn(x1, x2, . . . ,xn−1, xn)       

Marzouk et al. MIT 11 / 40

How to compute the sparsity pattern

3 2 5 1 4

G5

3 2 5 1 4

G4

3 2 5 1 4

G3

5 2 4 1 3

G2

◮ Compute marginal graphs: Gi−1 is obtained from Gi by removing

node i and by turning its neighborhood into a clique (like variable elimination)

◮ Sparsity of inverse transport: the i-th

component of S can depend, at most, on the variables in a neighborhood of node i in Gi

◮ Sparsity depends on the ordering of the variables

(similar heuristics sparse Cholesky)

Pkj = ∂xjSk

Marzouk et al. MIT 12 / 40

Decomposable transport maps

◮ Definition: a decomposable transport is a map T = T1 ◦ · · · ◦ Tk

that factorizes as the composition of finitely many maps of low effective dimension that are triangular (up to a permutation), e.g., T(x) =         A1(x1, x2, x3) B1(x2, x3) C1(x3) x4 x5 x6        

• T1
• 

       x1 A2(x2, x3, x4, x5) B2(x3, x4, x5) C2(x4, x5) D2(x5) x6        

• T2
• 

       x1 x2 x3 A3(x4) B3(x4, x5) C3(x4, x5, x6)        

• T3

◮ Theorem [SBM 2018]: Decomposable graphical models for π lead to

decomposable direct maps T, provided that η(x) =

i η(xi)

Marzouk et al. MIT 13 / 40

Decomposable transport maps

◮ Example graph decomposition V = (A, S, B) ◮ Effective dimension of each component map is |A ∪ S|

A B S

2 3 1 6 4 5

Marzouk et al. MIT 14 / 40

Transport maps and graphical models Key message

◮ Inverse maps: enforce sparsity in the approximation space S△, i.e.,

in solving minS∈S△ DKL( π || S−1

♯

η )

◮ Can also use for structure learning in non-Gaussian graphical models

[Morrison/Baptista/M NIPS 2017]

◮ Direct maps: enforce decomposable structure in the approximation

space T△, i.e., when solving minT∈T△ DKL( T♯η || π )

◮ A general tool for modeling and computation with non-Gaussian

Markov random fields

Marzouk et al. MIT 15 / 40

Transport maps and state-space models

Z0 Z1 Z2 Z3 ZN Y0 Y1 Y2 Y3 YN

Θ

◮ In many situations, elements of the composition

T = T1 ◦ T2 ◦ · · · ◦ Tk can be constructed sequentially

◮ Yields new algorithms for smoothing and and joint state-parameter

inference in state-space models [SBM 2018; Houssineau, Jasra, Singh

2018]

◮ Recent offshoot: nonlinear ensemble filtering [arXiv:1907.00389]

Marzouk et al. MIT 16 / 40

Remainder of the talk: two vignettes

1 Exploiting the low-dimensional structure of transport ◮ Sparsity ◮ Decomposability ◮ Low rank 2 What to do with approximate maps? ◮ Preconditioned/adaptive MCMC Marzouk et al. MIT 17 / 40

Vignette #1: low rank structure

◮ Let U = [Ur U⊥] ∈ Rn×n be a unitary matrix, with Ur ∈ Rn×r. A

lazy map T : Rn → Rn takes the form: T(z) = Urτ(z1, . . . , zr) + U⊥z⊥ for some diffeomorphism τ : Rr → Rr.

◮ Map T ∈ Tr(U) departs from the identity only on an r-dimensional

subspace

◮ Proposition: For any lazy map T ∈ Tr(U), there exists a strictly

positive function f : Rr → R+ such that T♯η(x) = f (U⊤

r x) η(x),

for all x ∈ Rn where η = N(0, In). Conversely, any density of the form f (U⊤

r x) η(x) for some f : Rr → R+ admits a lazy map

representation.

Marzouk et al. MIT 18 / 40

Low-rank structure

Why would such structure (approximately) appear?

◮ Bayesian inverse problems: data only partially informative; posterior

departs from the prior primarily on a low-dimensional subspace.

◮ Formalized by likelihood-informed subspace [Cui et al. 2014]; also,

active subspace [Constantine et al. 2015], and recent refinements/connections [Zahm et al. 2018].

Marzouk et al. MIT 19 / 40

Error bound and subspace

How to find a good Ur?

◮ Define

Hπ := ∇ log π η

• ∇ log π

η ⊤ dπ .

◮ Let (λi, ui) be the ith eigenpair of Hπ and put Ur = [u1 u2 · · · ur]. ◮ Theorem [Zahm et al. 2018]:

DKL(π||T ⋆

♯ η) ≤ 1

2(λr+1 + . . . + λd). where T ⋆

♯ η = f ⋆(U⊤ r x)η(x) and f ⋆(zr) = EX∼η

• π(X)

η(X)|U⊤ r X = zr

• .

◮ Good approximation when the spectrum of Hπ decays quickly ◮ Uses a ridge approximation of dπ/dη (e.g., the likelihood), with

• ptimal profile function f ⋆

Marzouk et al. MIT 20 / 40

Layers of lazy maps

◮ What if (λi) do not decay quickly? What if we are limited to small r? ◮ Answer: layers of lazy maps, via a greedy construction

◮ Given (π, η, r1): compute Hπ and construct a first lazy map T1 ◮ Pull back π by T1: π2 := (T −1

1 )♯π

◮ Given (π2, η, r2): compute Hπ2 and construct a next lazy map T2 . . . ◮ Generic iteration: at stage ℓ, build a lazy map to the pullback

πℓ := (T1 ◦ T2 ◦ · · · ◦ Tℓ−1)−1

♯ π

◮ Stop when 1

2 Tr(Hπℓ) < ǫ

Inputs Lazy map Outputs Rotation Rotation Lazy map Rotation

Marzouk et al. MIT 21 / 40

Layers of lazy maps

Example: rotated “banana” target distribution, r = 1 maps

Target π T♯

1π

T♯

2π

T♯

3π

T♯

5π

T♯

8π

Marzouk et al. MIT 22 / 40

Example: log-Gaussian Cox process

1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8

Field Λ Λ Λ⋆ and observations y⋆ Realizations of Λ Λ Λ ∼ πΛ

Λ Λ|y⋆

Marzouk et al. MIT 23 / 40

Example: log-Gaussian Cox process

◮ Parameter dimension n = 4096, 30 observations; fixed ranks r

1 2 3 4 5 6 7 Lazy iteration 100 101 102 103

1 2Tr(H )

Lazy rank: 1 Lazy rank: 3 Lazy rank: 5

Convergence

10 20 30 Eigenvalues 10

3

10

1

101 103 Eig(H0) Eig(H1) Eig(H2) Eig(H3) Eig(H4)

Spectrum of Hπℓ

Marzouk et al. MIT 24 / 40

Example: elliptic PDE Bayesian inverse problem

• ∇ · (eκ(x)∇u(x)) = 0,

for x ∈ D := [0, 1]2 , u(x) = 0 for x1 = 0, u(x) = 1 for x1 = 1,

∂u(x) ∂n

= 0 for x2 ∈ {0, 1}

◮ Infer κ(x), discretized with n = 2601 parameters; 81 observations;

lazy maps of r ≤ 4 and polynomial degree up to 2

1

u(x) and observations

2 4 6 8 10 Lazy iteration 100 101 102 103 104 105

1 2Tr(H ) 1 2V[log /T

]

Convergence

3 2 1 1 2

Posterior realizations

• f κ(x)

Marzouk et al. MIT 25 / 40

Vignette #2: preconditioning

What to do when T♯η = π?

◮ Maybe close enough? Can evaluate variance diagnostic

Varη[log(η/T −1

♯

¯ π)], bound Tr(HT −1

♯

π), etc. ◮ Enrich T, e.g., add a layer or expand T h △ in the given layer ◮ Work on the pullback: treat T −1 ♯

π with an asymptotically exact scheme, e.g., Markov chain Monte Carlo

Marzouk et al. MIT 26 / 40

Vignette #2: preconditioning

What to do when T♯η = π?

◮ Maybe close enough? Can evaluate variance diagnostic

Varη[log(η/T −1

♯

¯ π)], bound Tr(HT −1

♯

π), etc. ◮ Enrich T, e.g., add a layer or expand T h △ in the given layer ◮ Work on the pullback: treat T −1 ♯

π with an asymptotically exact scheme, e.g., Markov chain Monte Carlo One possible construction: transport-accelerated MCMC

◮ Transport map “preconditions” MCMC target; use MCMC iterates in

maps-from-samples construction

◮ Can be understood in the framework of adaptive MCMC

Marzouk et al. MIT 26 / 40

Preconditioning MCMC

◮ Effective MCMC proposal = adapted to the target

◮ Can we transform proposals or, equivalently, targets for better

sampling?

Marzouk et al. MIT 27 / 40

Recall maps-from-samples construction

◮ Suppose we have samples xi from target π

• S ∈ arg max

S∈Sh

△

1 M

M

• i=1

log S−1

♯

η (xi), η = N(0, In),

◮ Easy to solve: convex and separable, regardless of form of π ◮ View

S♯π as the “preconditioned” target

◮ In the MCMC setting, {xi}M

i=1 comprises dependent MCMC samples

◮

S♯π may be far from standard Gaussian for small M and/or crude Sh

△

Marzouk et al. MIT 28 / 40

• Ingredient #1: static map

– Perform MCMC in the reference space, on the “preconditioned” density – Simple proposal in reference space (e.g., random walk) corresponds to a more complex/tailored proposal on target

r θ ˜ p(r) ˜ T(θ) π(θ)

Map-accelerated MCMC

S(x)

π(x)

x

S]π

• Ingredient #1: static map

– Perform MCMC in the reference space, on the “preconditioned” density – Simple proposal in reference space (e.g., random walk) corresponds to a more complex/tailored proposal on target

α = π(S −1( ′ r ))| ∇S −1 | ′

r qr(r | ′

r ) π(S −1(r))| ∇S −1 |r qr( ′ r |r)

simple proposal qr on pushforward of target through map

Map-accelerated MCMC

r θ ˜ p(r) S(θ) π(θ) θ r r0 qr(r0|r) θ0

S]π

S(x)

π(x)

x

x0

x

• Ingredient #1: static map

– Perform MCMC in the reference space, on the “preconditioned” density – Simple proposal in reference space (e.g., random walk) corresponds to a more complex/tailored proposal on target

α = π(S −1( ′ r ))| ∇S −1 | ′

r qr(r | ′

r ) π(S −1(r))| ∇S −1 |r qr( ′ r |r)

more complex proposal, directly on target distribution

Map-accelerated MCMC

r θ ˜ p(r) S(θ) π(θ) θ r r0 qr(r0|r) θ0

S]π

S(x)

π(x)

x

x0

x

• Ingredient #2: adaptive map

– Update the map with each MCMC iteration: more samples, more accurate , better – Adaptive MCMC [Haario 2001, Andrieu 2006], but with nonlinear transformation to capture non-Gaussian structure S

r θ ˜ Tk(θ) ˜ pk(r) π(θ)

Map-accelerated MCMC

Sk(θ)

Eπ

π(x)

x

Sk(x)

(Sk)]π

r θ ˜ Tk+1(θ) π(θ) ˜ pk+1(r)

• Ingredient #2: adaptive map

– Update the map with each MCMC iteration: more samples, more accurate , better – Adaptive MCMC [Haario 2001, Andrieu 2006], but with nonlinear transformation to capture non-Gaussian structure

Map-accelerated MCMC

π(x)

x

(Sk+1)]π

Sk+1(x)

S Eπ

• Ingredient #3: global proposals

– If the map becomes sufficiently accurate, would like to avoid random-walk behavior

Map-accelerated MCMC

reference RW proposal mapped RW proposal

• Ingredient #3: global proposals

– If the map becomes sufficiently accurate, would like to avoid random-walk behavior

Map-accelerated MCMC

reference independence proposal mapped independence proposal

• Ingredient #3: global proposals

– If the map becomes sufficiently accurate, would like to avoid random-walk behavior – Solution: delayed rejection MCMC [Mira 2001] – First proposal = independent sample from η (global, more efficient); second proposal = random walk (local, more robust)

• Entire scheme is provably ergodic with respect to the exact

posterior measure [Parno & M, SIAM JUQ 2018]

– Requires enforcing some regularity conditions on maps, to preserve tail behavior of transformed target

Map-accelerated MCMC

Example: biological oxygen demand model

◮ Likelihood model:

d = θ1(1 − exp (−θ2x)) + ǫ ǫ ∼ N

• 0, 2 × 10−4

◮ 20 noisy observations at

x = 5 5, 6 5, . . . , 25 5

• ◮ Degree-three polynomial map

True posterior density

Marzouk et al. MIT 29 / 40

Transformed distribution

Original posterior π Pushforward of posterior through learned map, S♯π

Marzouk et al. MIT 30 / 40

Results: ESS per computational effort

T M + D R T M + N U T S T M + L A N N U T S D R A M 100 200 300 161 57 179 1.4 5.8 ESS/(1,000 Evaluations) – θ1 T M + D R G T M + D R L T M + M I X N U T S D R A M 500 1,000 1,500 1468 487 1495 57 127 ESS/second– θ1

Marzouk et al. MIT 31 / 40

Example #2: predator-prey model

◮ Six parameter ODE population model

dP dt = rP

• 1 − P

K

• − s PQ

a + P dQ dt = u PQ a + P − vQ

◮ Five noisy observations of both populations ◮ Infer 6 parameters + 2 initial values; uniform priors

Marzouk et al. MIT 32 / 40

Predator-prey model: chains

1 2 3 P(0) DRAM sMMALA AMALA 1 2 3 MCMC Step P(0) TMDRG MCMC Step TMDRL MCMC Step TMRWM

Marzouk et al. MIT 33 / 40

Example: maple sap dynamics model

◮ Coupled PDE system for

ice, water, and gas locations [Ceseri & Stockie 2013]

◮ Measure gas pressure in

vessel

◮ Infer 10 physical model

parameters

◮ Very challenging

posterior!

R +R+r R s

x

R+R−r s R+2R water

fiber/ vessel vessel wall

f

fiber

R

v

sgi siw gas ice water gas

gi iw f f v f v f v

r R

Image from Ceseri and Stockie, 2013

Marzouk et al. MIT 34 / 40

Maple posterior distribution

θ1 θ2 θ3 θ4 θ5 θ6 θ7 θ8 θ9 θ10

−1 −0.5 0.5 1 2 θ1 θ7 −1 −0.5 0 0.5 1 −0.6 −0.4 −0.2 0.2 0.4 θ3 θ6

Marzouk et al. MIT 35 / 40

Results: ESS per computational effort

T M + D R G T M + D R L T M + M I X D R A M 2 4 6 8 10 5.7 10 2.9 0.6 ESS/(10,000 Evaluations) T M + D R G T M + D R L T M + M I X D R A M 5 10 15 20 25 18 26 7.1 2.3 ESS/(1000 seconds)

Marzouk et al. MIT 36 / 40

Comments on MCMC with transport maps

Useful characteristics of the algorithm:

◮ Map construction is easily parallelizable ◮ Requires no gradients from posterior density

Generalizes many current MCMC techniques:

◮ Adaptive Metropolis: map enables non-Gaussian proposals and a

natural mixing between local and global moves

◮ Manifold MCMC [Girolami & Calderhead 2011]: map also defines a

Riemannian metric

Marzouk et al. MIT 37 / 40

Looking to higher dimensions: regularized estimation of S

For simplicity, consider map components Sk(x) =

j βjψj(x1:k−1) + αkxk

• Sk ∈ arg

min

Sk∈Sh

△,k

1 N

N

• i=1

1 2Sk(xi)2 − log ∂kSk(xi)

• + λNβ1

Assume sub-Gaussian π and basis functions ψj(x)

Theorem [BZM]

For polynomial maps of degree m with sparsity s, with high probability Eπ

• DKL
• π(xk|x1:k−1) ||

S♯

kη

• s2m log k

N Takeaways

◮ Accurate estimation is feasible in high dimensions with N ≪ k ◮ From factorization property of density, error in conditionals ensures

DKL(π|| S♯η) d

• s2m log d

N

Marzouk et al. MIT 38 / 40

Conclusions

◮ Central idea: characterize complex/intractable distributions by

constructing deterministic couplings

◮ Many kinds of low-dimensional structure (non-exhaustive):

◮ Sparse maps, decomposable maps ◮ Low rank structure (lazy maps)

◮ Exploiting the pullback distribution

◮ Compositions of approximate maps, constructed greedily ◮ Use approximate maps to precondition other sampling or cubature

schemes

Marzouk et al. MIT 39 / 40

Conclusions

Extensions and ongoing work:

◮ Using sparse grids or QMC for map construction ◮ Preconditioning other sampling schemes: computational tradeoffs,

impact on target geometry and regularity

◮ Zoo of map parameterizations and their approximation properties ◮ Additional varieties of low-dimensional structure ◮ Sequential variational inference with maps: filtering, smoothing,

parameter estimation

◮ Maps from samples:

◮ Nonparametric estimation schemes ◮ As a tool for ensemble filtering [arXiv:1907.00389] or approximate

Bayesian computation. . .

Marzouk et al. MIT 39 / 40

Conclusions

Extensions and ongoing work:

◮ Using sparse grids or QMC for map construction ◮ Preconditioning other sampling schemes: computational tradeoffs,

impact on target geometry and regularity

◮ Zoo of map parameterizations and their approximation properties ◮ Additional varieties of low-dimensional structure ◮ Sequential variational inference with maps: filtering, smoothing,

parameter estimation

◮ Maps from samples:

◮ Nonparametric estimation schemes ◮ As a tool for ensemble filtering [arXiv:1907.00389] or approximate

Bayesian computation. . .

Thanks for your attention!

Marzouk et al. MIT 39 / 40

References

◮ A. Spantini, R. Baptista, Y. Marzouk. “Coupling techniques for nonlinear ensemble

filtering.” arXiv:1907.00389.

◮ D. Bigoni, O. Zahm, A. Spantini, Y. Marzouk. “Greedy inference with layers of

lazy maps.” arXiv:1906.00031.

◮ O. Zahm, T. Cui, K. Law, A. Spantini, Y. Marzouk. “Certified dimension

reduction in nonlinear Bayesian inverse problems.” arXiv:1807.03712.

◮ A. Spantini, D. Bigoni, Y. Marzouk. “Inference via low-dimensional couplings.”

JMLR 19(66): 1–71, 2018.

◮ M. Parno, Y. Marzouk, “Transport map accelerated Markov chain Monte Carlo.”

SIAM JUQ 6: 645–682, 2018.

◮ G. Detomasso, T. Cui, A. Spantini, Y. Marzouk, R. Scheichl, “A Stein variational

Newton method.” NeurIPS 2018.

◮ R. Morrison, R. Baptista, Y. Marzouk. “Beyond normality: learning sparse

probabilistic graphical models in the non-Gaussian setting.” NeurIPS 2017.

◮ Y. Marzouk, T. Moselhy, M. Parno, A. Spantini, “An introduction to sampling via

measure transport.” Handbook of Uncertainty Quantification, R. Ghanem, D. Higdon, H. Owhadi, eds. Springer (2016). arXiv:1602.05023.

◮ General python code at http://transportmaps.mit.edu

Marzouk et al. MIT 40 / 40