Stratification for Rare Event Problems

Department of Statistics & 

James Franck Institute 

University of Chicago

April 27th, 2017

Jonathan Weare 

Aaron Dinner, Andrey Kravstov, Charles Matthews, Jonathan Mattingly, 

Jeremy Tempkin, Erik Thiede, and Brian Van Koten
The rare event problem

A B

transition 

barrier

q

Energy

G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187–199 (1977) 

C. Pangali, M. Rao, and B. J. Berne, J. Chem. Phys. 71, 2975-2981 (1979) 

E. Thiede, B. Van Koten, A. Dinner, and JW, J. Chem. Phys. 145(8), 084115 (2016)

rV (Xt) dt + p2" dWt

V

Can sample 

V

⇡ / eV /"

but very slow convergence 

(exponentially in ) for most averages

dXt = rV (Xt) dt + p2" dWt

V

⇡ / eV /"

"1

Y (t+1) = X(t) + N (0, )

X(t+1) =

(

Y (t+1), with probability 1 ^ e(V (Y (t+1))

X(t), otherwise

 ̃ij (m)

J(t) 2 {0, 1, 2,...,n}

eF= X

n

zi heW(t)ii

Y (t+1) = X(t) + N (0, )

X(t+1) =

(

Y (t+1), w.p. 1 ^ ⇡(Y (t+1))

⇡(X(t))

X(t), otherwise

 ̃ij (m)

J(t) 2 {0, 1, 2,...,n}

eF= X

n

i=1

zi heW(t)ii

F

E

h

W(t)

i

with e.g. Metropolis:
Umbrella Sampling

G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187–199 (1977) 

C. Pangali, M. Rao, and B. J. Berne, J. Chem. Phys. 71, 2975-2981 (1979) 

E. Thiede, B. Van Koten, A. Dinner, and JW, J. Chem. Phys. 145(8), 084115 (2016)

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

 i  0

Z

f(x)⇡(dx) = X

i

zi hfii

E

"

X

t=0

f(t, X(t))

#

= X

i

z

J = j

X

i

 i = 1

Z

f(x)⇡(dx) = X

i

zi hf

1. Introduce “bias” functions

Target

distribution

Target distribution

x

x

0

-2

-4

-6

-8

2

Biasing

0

-2

-4

-6

-8

2

ln(π)

ln(π)

ln(π1)

ln(π2)

ln(π3)

ln(π4)

ln(π5)

ln(ψ5)

ln(ψ4)

ln(ψ3)

ln(ψ2)

ln(ψ1)
Umbrella Sampling

G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187–199 (1977) 

C. Pangali, M. Rao, and B. J. Berne, J. Chem. Phys. 71, 2975-2981 (1979) 

E. Thiede, B. Van Koten, A. Dinner, and JW, J. Chem. Phys. 145(8), 084115 (2016)

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

 i  0

Z

f(x)⇡(dx) = X

i

zi hfii

E

"

X

t=0

f(t, X(t))

#

= X

i

z

J = j

X

i

 i = 1

Z

f(x)⇡(dx) = X

i

zi hf

1. Introduce “bias” functions

2. Define “restrained” distributions

{zi}

zj = X

i

Fij zi, X

i

zi = 1

Fij =

Z

 j (x) ⇡i(dx)

⇡i / i ⇡

 log i

Z

f(x)⇡(dx) = X

i

zi hfii

Target

distribution

Target distribution

x

x

0

-2

-4

-6

-8

2

Biasing

0

-2

-4

-6

-8

2

ln(π)

ln(π)

ln(π1)

ln(π2)

ln(π3)

ln(π4)

ln(π5)

ln(ψ5)

ln(ψ4)

ln(ψ3)

ln(ψ2)

ln(ψ1)
Umbrella Sampling

G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187–199 (1977) 

C. Pangali, M. Rao, and B. J. Berne, J. Chem. Phys. 71, 2975-2981 (1979) 

E. Thiede, B. Van Koten, A. Dinner, and JW, J. Chem. Phys. 145(8), 084115 (2016)

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

 i  0

Z

f(x)⇡(dx) = X

i

zi hfii

E

"

X

t=0

f(t, X(t))

#

= X

i

z

J = j

X

i

 i = 1

Z

f(x)⇡(dx) = X

i

zi hf

1. Introduce “bias” functions

2. Define “restrained” distributions

{zi}

zj = X

i

Fij zi, X

i

zi = 1

Fij =

Z

 j (x) ⇡i(dx)

⇡i / i ⇡

 log i

Z

f(x)⇡(dx) = X

i

zi hfii

3. Define normalization constants

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

zi =

Z

 i(x) ⇡(dx)

Z

f(x)⇡(dx) = X

i

zi hfii

Target

distribution

Target distribution

x

x

0

-2

-4

-6

-8

2

Biasing

0

-2

-4

-6

-8

2

ln(π)

ln(π)

ln(π1)

ln(π2)

ln(π3)

ln(π4)

ln(π5)

ln(ψ5)

ln(ψ4)

ln(ψ3)

ln(ψ2)

ln(ψ1)
4. Assemble general averages

eF= X

n

i=1

zi heW(t)ii

eF= E

h

eW(t)

i

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

Z

f(x) ⇡(dx) = X

i

zi

Z

f(x) ⇡i(dx)

Z

f(x)⇡(dx) = X

i

zi hfii

Umbrella Sampling

G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187–199 (1977) 

C. Pangali, M. Rao, and B. J. Berne, J. Chem. Phys. 71, 2975-2981 (1979) 

E. Thiede, B. Van Koten, A. Dinner, and JW, J. Chem. Phys. 145(8), 084115 (2016)

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

 i  0

Z

f(x)⇡(dx) = X

i

zi hfii

E

"

X

t=0

f(t, X(t))

#

= X

i

z

J = j

X

i

 i = 1

Z

f(x)⇡(dx) = X

i

zi hf

1. Introduce “bias” functions

2. Define “restrained” distributions

{zi}

zj = X

i

Fij zi, X

i

zi = 1

Fij =

Z

 j (x) ⇡i(dx)

⇡i / i ⇡

 log i

Z

f(x)⇡(dx) = X

i

zi hfii

3. Define normalization constants

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

zi =

Z

 i(x) ⇡(dx)

Z

f(x)⇡(dx) = X

i

zi hfii

Target

distribution

Target distribution

x

x

0

-2

-4

-6

-8

2

Biasing

0

-2

-4

-6

-8

2

ln(π)

ln(π)

ln(π1)

ln(π2)

ln(π3)

ln(π4)

ln(π5)

ln(ψ5)

ln(ψ4)

ln(ψ3)

ln(ψ2)

ln(ψ1)
Umbrella Sampling

i.e. the weights solve the eigenproblem

V

⇡ / eV /"

"1

{zi}

zj = X

i

Fij zi, X

i

zi = 1

Fij =

with

V

⇡ / eV /"

"1

{zi}

zj = X

i

Fij zi, X

i

zi = 1

Fij =

Z

 j (x) ⇡i(dx)

⇡/ ⇡

E. Thiede, B. Van Koten, A. Dinner, and JW, J. Chem. Phys. 145(8), 084115 (2016)

⇡ / eV /"

"1

{zi}

zT= zTF, X

i

zi = 1

Fij =

Z

 j (x) ⇡i(dx)

X

 i = 1

"

X

t=0

#

X

i

J = j

Z

f(x) ⇡(dx) = X

i

zi

Z

f(x) ⇡i(dx)

Z

f(x)⇡(dx) = X

i

zi hfii

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

F

zj =

Z

 j (x)⇡(dx) = X

i

zi

Z

 j (x)⇡i(dx)

2

so 
1. Use you favorite MC scheme to sample from 

2. Estimate the entries of and each 

3. Solve the eigenproblem for 

4. Assemble 

Z

f(x)⇡(dx) = X

i

z

⇡i

zj =

Z

 j (x)⇡(dx) = X

i

zi

Z

E

"

X

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

F

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

Zh

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

Z

f(x)⇡i(dx)

zj =

Z

 j (x)⇡(dx) = X

i

zi

Z

 j 

Z

f(x)⇡(dx) = X

i

zi hfii

Z

f(x)⇡i(dx)

z

zj =

Z

 j (x)⇡(dx) = X

i

zi

Z

 j (x)⇡i(d

hi

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

Z

f(x) ⇡(dx) = X

i

zi

Z

f(x) ⇡i(dx)

Z

f(x)⇡(dx) = X

i

zi hfii

Umbrella Sampling

E. Thiede, B. Van Koten, A. Dinner, and JW, J. Chem. Phys. 145(8), 084115 (2016)
A detailed analysis

We establish a central limit theorem with a detailed 

expression for the asymptotic variance along with a 

general bound:

E. Thiede, B. Van Koten, A. Dinner, and JW, J. Chem. Phys. 145(8), 084115 (2016)

F

= X

n

i=1

zi heW(t)ii

eF= E

h

eW(t)

i

hfii

zT= zTG + aT

⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

Pi [tj < ti]

i

is the probability that the chain defined by and starting from 

index and hitting index before returning to 

eF= X

n

i=1

zi heW(t)ii

eF= E

h

eW(t)

i

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

F

⌦

{}

eF= X

n

i=1

zi heW(t)ii

eF= E

h

eW(t)

i

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

i

j

J(t) 2 {0, 1, 2,...,n}

eF= X

n

i=1

zi heW(t)ii

eF= E

h

eW(t)

i

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

i

j

i

{}

eF= X

n

i=1

zi heW(t)i

eF= E

h

eW(t)

i

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi 

⇡i / (⇡)

↵i

i

j

contribution from 

restricted average of 

hi

hfii

zT = zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

f

j

i

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

contribution from 

estimate of 

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

F

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

2(f)  2

X

i

1

i

8

>><

>>:

⌧iz2

i + ⇡[|f|]

2tr(Ri

)

X

j6=i

Fij 6=0

var⇡i ( j )

Pi[tj < ti]

9

>>=

>>;

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k)

f

i

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

2(f)  2

X

i

1

i

8

>><

>>:

⌧iz2

i + ⇡[|f|]

2tr(Ri

)

X

j6=i

Fij 6=0

var⇡i ( j )

Pi[tj < ti]

9

>>=

>>;

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k)

f

i

⌦

is the asymptotic covariance of the estimators of and

= zTG + aT

X(t))

#

= X

i

zi hfii

i / (⇡)

↵i

+ ⇡[|f|]

2tr(Ri

)

X

j6=i

Fij 6=0

var⇡i ( j )

Pi[tj < ti]

9

>>=

>>;

i

jk

r⇡i ( j )

pvar⇡i ( k)

f

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

2(f)  2

X

i

1

i

8

>><

>>:

⌧iz2

i + ⇡[|f|]

2tr(Ri

)

X

j6=i

Fij 6=0

var⇡i ( 

Pi[tj < 

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k)

Fij

Fik

"

t=0

#

i

⇡i / (⇡)

↵i

2(f)  2

X

i

1

i

8

>><

>>:

⌧iz2

i + ⇡[|f|]

2tr(Ri

)

X

j6=i

Fij 6=

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k

Fij

Fik

i
sensitive to some entries of the transition matrix, but 

example arose from a problem in computational statis

the example to compare our results with previous work, 

bounds on absolute error summarized in [3].

6.1. Transition matrix and physical interpretatio

ample is a simple analogue of the dynamics of a single 

heat bath. Define V : R ! R by

V (x) = 1

4⇡

cos(4⇡x)

Take L 2 N, and let ⌦ := {1, 2,...,L} with periodic bo

take ⌦ := Z/LZ. Given V , we define a probability distr

⇡(i) := exp(LV (i/L))

PL

k=1 exp(LV (k/L))

The measure ⇡ is in detailed balance with the Markov 

matrix F 2 R⌦⇥⌦ defined by

Fii :=

1

2

✓ ⇡(i)

⇡(i  1) + ⇡(i) +

⇡(i)

⇡(i + 1) + ⇡(i)

◆

Fi,i+1 :=

1

2

⇡(i + 1)

⇡(i + 1) + ⇡(i) 

1(i 1)

ample is a simple analogue of the dynamics of a single particle in conta

heat bath. Define V : R ! R by

V (x) = 1

4⇡

cos(4⇡x)

Take L 2 N, and let ⌦ := {1, 2,...,L} with periodic boundary condition

take ⌦ := Z/LZ. Given V , we define a probability distribution on ⌦ by

⇡(i) := exp(LV (i/L))

PL

k=1 exp(LV (k/L)).

The measure ⇡ is in detailed balance with the Markov chain X having t

matrix F 2 R⌦⇥⌦ defined by

Fii :=

1

2

✓ ⇡(i)

⇡(i  1) + ⇡(i) +

⇡(i)

⇡(i + 1) + ⇡(i)

◆

for all i 2 ⌦,

Fi,i+1 :=

1

2

⇡(i + 1)

⇡(i + 1) + ⇡(i) for all i 2 ⌦,

Fi,i1 :=

1

2

⇡(i  1)

⇡(i  1) + ⇡(i) for all i 2 ⌦, and

Fij := 0 otherwise.

on the periodic 1D lattice

Consider the probability vector:

with 

6. The hilly landscape example

n this section, we discuss an example in which the invariant distribution is very

itive to some entries of the transition matrix, but insensitive to others. The

mple arose from a problem in computational statistical physics. We will use

example to compare our results with previous work, especially [2, 7, 13] and the

nds on absolute error summarized in [3].

Transition matrix and physical interpretation. Our hilly landscape ex

ple is a simple analogue of the dynamics of a single particle in contact with a

t bath. Define V : R ! R by

V (x) = 1

4⇡

cos(4⇡x)

e L 2 N, and let ⌦ := {1, 2,...,L} with periodic boundary conditions; that is,

⌦ := Z/LZ. Given V , we define a probability distribution on ⌦ by

⇡(i) := exp(LV (i/L))

PL

k=1 exp(LV (k/L)).

measure ⇡ is in detailed balance with the Markov chain X having transition

rix F 2 R⌦⇥⌦ defined by

Fii :=

1

2

✓ ⇡(i)

⇡(i  1) + ⇡(i) +

⇡(i)

⇡(i + 1) + ⇡(i)

◆

for all i 2 ⌦,

Fi,i+1 :=

1

2

⇡(i + 1)

⇡(i + 1) + ⇡(i) for all i 2 ⌦,

Fi,i1 :=

1

2

⇡(i  1)

⇡(i  1) + ⇡(i) for all i 2 ⌦, and

Fij := 0 otherwise

An illustrative example

A. Dinner, E. Thiede, B. Van Koten, and JW, SIMAX 36(3), 917-941 (2015)

Nearest neighbor Metropolis MCMC generates a chain with 

exponentially small spectral gap in 

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

2(f)  2

X

i

1

i

8

>><

>>:

⌧iz2

i + ⇡[|f|]

2tr(Ri

)

X

j6=i

Fij 6=0

var⇡i 

Pi[tj 

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k)

L

F

⇡
An illustrative example

A. Dinner, E. Thiede, B. Van Koten, and JW, SIMAX 36(3), 917-941 (2015)

Fij =

8

>>>><

>>>>:

1

2

⇡(i+1)

⇡(i)+⇡(i+1) , j = i + 1

1

2

⇡(i)

⇡(i)+⇡(i+1) , j = i  1

1

2 , j = i

0, otherwise

⌧

Fik

tr(Ri) ⇠

1

# of MCMC steps

Pi[tj < ti]  Fij

Consider umbrella sampling with

Then

and 

Fij =

8

>>>><

>>>>:

1

2

⇡(i+1)

⇡(i)+⇡(i+1) , j = i + 1

1

2

⇡(i)

⇡(i)+⇡(i+1) , j = i  1

1

2 , j = i

0, otherwise

zi =

1

2

(⇡(i) + ⇡(i + 1))

Fik

tr(Ri) ⇠

1

# of MCMC steps

:

zi =

1

2

(⇡(i) + ⇡(i + 1))

 i(j) = 1

2

(i(j) + i+1(j))

tr(Ri) ⇠

1

# of MCMC steps

Pi[tj < ti]  Fij

Fi,i±1 =

1

2

1

1 + exp ⇥LV  i±1

L



 V  i

L

⇤ 

1

2

1

1 + e

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

Fij =

8

>>>><

>>>>:

1

2

⇡(i+1)

⇡(i)+⇡(i+1) , j = i + 1

1

2

⇡(i)

⇡(i)+⇡(i+1) , j = i  1

1

2 , j = i

0, otherwise

zi =

1

2

(⇡(i) + ⇡(i + 1))

 i(j) = 1

2

(i(j) + i+1(j))

tr(Ri) ⇠

1

# of MCMC steps

⇡i(j) = ⇡(i)i(j) + ⇡(i + 1)i+1(j)

⇡(i) + ⇡(i + 1)

Pi[tj < ti]  Fij

Fi,i±1 =

1

2

1

1 

⇥

LV  i±1



V  i

⇤ 

1

2

1

1 + 

We can sample e.g. by Metropolis

with simple random walk proposals
The inequality:

and also that

An illustrative example

A. Dinner, E. Thiede, B. Van Koten, and JW, SIMAX 36(3), 917-941 (2015)

>>>>:

2 , j 

i

0, otherwise

zi =

1

2

(⇡(i) + ⇡(i + 1))









log ⇡(i)

⇡(i + 1)









= L









V

✓i + 1

L

◆

 V

✓ i

L

◆







 1

 i(j) = 1

2

(i(j) + i+1(j))

tr(Ri) ⇠

1

# of MCMC steps

⇡i(j) = ⇡(i)i(j) + ⇡(i + 1)i+1(j)

⇡(i) + ⇡(i + 1)

Pi[tj < ti]  Fij

Fi,i±1 =

1

2

1

1 + exp ⇥LV  i±1

L



 V  i

L

⇤ 

1

2

1

1 + e

⌦

implies that sampling of is not hard, i.e.:

p

⇡i ( j )

p

⇡i ( k)

L

Fik

⇡i

1 =

1

2

1

1 + exp ⇥LV  i±1

L



 V  i

L

⇤ 

1

2

1

1 + e

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

>:

0, otherwise

zi =

1

2

(⇡(i) + ⇡(i + 1))









log ⇡(i)

⇡(i + 1)









= L









V

✓i + 1

L

◆

 V

✓ i

L

◆







 1

 i(j) = 1

2

(i(j) + i+1(j))

tr(Ri) ⇠

1

# of MCMC steps in window i

⇡i(j) = ⇡(i)i(j) + ⇡(i + 1)i+1(j)

⇡(i) + ⇡(i + 1)

Pi[tj < ti]  Fij

Fi,i±1 =

1

2

1

1 + exp ⇥LV  i±1

L



 V  i

L

⇤ 

1

2

1

1 + e

zi =

1

2

(⇡(i) + ⇡(i + 1))









log ⇡(i)

⇡(i + 1)









= L









V

✓i + 1

L

◆

 V

✓ i

L

◆







 1

 i(j) = 1

2

(i(j) + i+1(j))

tr(Ri) ⇠

1

# of MCMC steps in window i

⇡i(j) = ⇡(i)i(j) + ⇡(i + 1)i+1(j)

⇡(i) + ⇡(i + 1)

Pi[tj < ti]  Fij

Fi,i±1 =

1

2

1

1 + exp ⇥±LV  i+1

L



 V  i

L

⇤ 

1

2

1

1 + e

⌦

Z
Fik

tr(Ri) ⇠

1

# of MCMC steps

Pi[tj < ti]  Fij

=

1

2

1

1 + exp ⇥LV  i±1

L



 V  i

L

⇤ 

1

2

1

1 + e

⌦

1

And implies that 

So

⇡i / (⇡)

↵i

2(f)  2

X

i

1

i

8

>><

>>:

⌧iz2

i + ⇡[|f|]

2tr(Ri

)

X

j6=i

Fij 6=0

var⇡i ( j )

Pi[tj < ti]

9

>>=

>>;

X

j6=i

Fij 6=0

var⇡i ( j )

Pi[tj < ti]

⇠ 1

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k)

L

Fik

tr(Ri) ⇠

1

# of MCMC steps

"

X

t=0

f(, )

#

X

i

hfi

⇡i / (⇡)

↵i

2(f)  2

X

i

1

i

8

>><

>>:

⌧iz2

i + ⇡[|f|]

2tr(Ri

)

X

j6=i

Fij 6=0

var⇡i ( j )

Pi[tj < ti]

9

>>=

>>;

X

j6=i

Fij 6=0

var⇡i ( j )

Pi[tj < ti]

⇠ 1

2(f) ⇠

L

# of MCMC steps per window

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k)

L

Fik

and total cost to achieve fixed accuracy scales about like 

2(f)  2

X

i

1

i

><

>>:

⌧iz2

i + ⇡[|f|]

2tr(Ri

)

X

j6=i

Fij 6=

X

j6=i

Fij 6=0

var⇡i ( j )

Pi[tj < ti]

⇠ 1

2(f) ⇠

L

# of MCMC steps per w

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k

L2

Fik

t(Ri) 

1

i.e. not exponentially!

An illustrative example

A. Dinner, E. Thiede, B. Van Koten, and JW, SIMAX 36(3), 917-941 (2015)
US can reduce cost to sample from 

exponential in to algebraic.

V

⇡ / eV /"

⇡ / eV /"

"1

A detailed analysis

The restrained distributions can be 

much easier to sample. E.g. if 

is constant on a strip of width 

and if is sampled using reflected 

overdamped Langevin, then 

But: You have to think about how 

the eigenproblem amplifies errors 

in 

E. Thiede, B. Van Koten, A. Dinner, and JW, in preparation

tr(Ri) ⇠

1

# of MCMC steps in windo

⇡i(j) = ⇡(i)i(j) + ⇡(i + 1)i+1(j

⇡(i) + ⇡(i + 1)

Pi[tj < ti]  Fij

Fi,i±1 =

1

2

1

1 + exp ⇥±LV  i+1

L



 V  i

L

⇤ 

Zh

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

Zh

⌦

Z

f(x)⇡(dx) = X

i

 i

i,i±1 

2

1 + exp ⇥±LV  i+1

L



Zh

⌦

Z

f(x)⇡(dx) = X

i

h

Fi,i±1 

2

1 + exp ⇥±LV  i+1

L



 V  i

L

⇤ 

21 + 

Zh

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

⇡i

⇡i(j) = ⇡(i)i(j) + ⇡(i + 1)i+1(j)

⇡(i) + ⇡(i + 1)

Pi[tj < ti]  Fij

Fi,i±1 =

1

2

1

1 + exp ⇥±LV  i+1

L



 V  i

L

⇤ 

1

2

1

1 + e

Zh

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

F

Y (t+1) = X(t) + N (0, )

X(t+1) =

(

Y (t+1), w.p. 1 ^ ⇡(Y (t+1))

⇡(X(t))

X(t), otherwise

2

h(f) 

⇤h2varh(f)

"

exp ✓maxZh V  minZh V

"

◆

 ̃ij (m)

J(t) 2 {0, 1, 2,...,n}

eF= X

n

i=1

zi heW(t)ii

eF= E

h

eW(t)

i
A detailed analysis

US can reduce the cost to compute tail 

probabilities from exponential in 

to algebraic.

"1

P [X  M]

ai = P[J(0) = i]

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

"1

M

ai = P[J(0) = i

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

J = j

Instead of trying to estimate the 

small probability with relative 

accuracy, US estimates the 

relatively large entries in . 

But: Again, you have to think 

about how the eigenproblem 

amplifies errors in .

eF= X

n

i=1

zi heW(t)ii

eF= E

h

eW(t)

i

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

F

⌦

Z

f(x)⇡(dx) = Xzhfi

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi 

⇡i / (⇡)

↵i

F

⌦

Z

f(x)⇡(dx) = X

i

zi hfi

E. Thiede, B. Van Koten, A. Dinner, and JW, in preparation
A key ingredient

We need to understand the error in the 

resulting from sampling error in a stochastic matrix 

⇡ / e

/

"1

{zi}

zT= zTF, X

i

zi = 1

Fij =

Z

 j (x) ⇡i(dx)

X

i

 i = 1

 log i

⇡i / (⇡)

↵i

zi

⌦

Z

f(x)⇡(dx) = X

i

zi 

"

X

t=0

#

⇡i / (⇡)

↵

F

⌦

Z

f(x)⇡(dx) = X

Existing perturbation bounds for Markov matrices 

blow up as (one over) the spectral gap of 

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

⇡i / (⇡)

↵i

F

⌦

Z

f(x)⇡(dx) = Xzi 

In US applications the spectral gap of 

is typically extremely small. 

A. Dinner, E. Thiede, B. Van Koten, and JW, SIMAX 36(3), 917-941 (2015)
A key ingredient

0 10 20 30 40 50 60 70 80

L

0

2

4

6

8

10

12

14

16

18

log

i(F)

log 3(F)

log 7(F)

We are far from the first to consider this problem

Existing bounds blow up as the spectral gap decreases

A. Dinner, E. Thiede, B. Van Koten, and JW, SIMAX 36(3), 917-941 (2015)

Some bounds of the form (from Chao and Meyer 2001):

X=

(

()

X(t), otherwise

2

h(f) 

⇤h2varh(f)

"

exp ✓maxZh V  minZh V

"

◆

 ̃ij (m)

kz ̃  zk  ikF ̃  Fk

J(t) 2 {0, 1, 2,...,n}

eF= X

n

i=1

zi heW(t)ii

eF= E

h

eW(t)

i

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii
A key ingredient

We establish new more detailed bounds that show 

that error can be small even when spectral gap is small 

e e.g. [?]). Very generally, existing results tell us that, when the mixing of 

low, the invariant measure will be highly sensitive to general perturbations 

the process. Unfortunately, existing results are not suciently detailed to 

rturbations that do not lead to large errors and those that do. In many applic

are crucial. We have recently introduced very detailed perturbation bound

ov processes which revealed that, for example, that perturbations resulting fro

a widely used family of Monte Carlo methods (umbrella sampling) do not typ

tive errors even when the spectral gap of the relevant matrix is tiny [?, ?]. Mo

at, if z and  ̃z are the invariant probability vectors for stochastic matrices F 

substochastic matrix S, the entries of F and F ̃ are bounded below by the co

then

max

k |log zk  log  ̃zk||  X

i6=j

|log(Fij + ij )  log(F ̃ij + ij )|

Pi(⌧j < min{⌧i, ⌧})  Sij where ⌧i is the return time to state i for the absorb

ned by S, and ⌧ is the time to absorption for that chain. The first importan

ity is that it bounds a form of relative error in the entries of the invariant 

s is crucial in the rare event context where the goal is to compute small p

xample that removing the ij (which are all non-negative) from the right hand 

hows that the relative error in the invariant probability vector is bounded by 

e non-zero entries of F. It is not dicult to construct simple stochastic matric

itlllbt thtiiF ld

of a Markov process resulting from misspecification of the p

problem (see e.g. [?]). Very generally, existing results tell u

process is slow, the invariant measure will be highly sensit

rameters of the process. Unfortunately, existing results are 

between perturbations that do not lead to large errors and t

distinctions are crucial. We have recently introduced very 

state Markov processes which revealed that, for example, tha

errors from a widely used family of Monte Carlo methods (um

in large relative errors even when the spectral gap of the relev

we prove that, if z and  ̃z are the invariant probability vecto

if, for some substochastic matrix S, the entries of F and F ̃ a

entry in S, then

max

k |log zk  log  ̃zk||  X

i6=j

|log(Fij + i

where ij = Pi(⌧j < min{⌧i, ⌧})  Sij where ⌧i is the return t

chain governed by S, and ⌧ is the time to absorption for th

this inequality is that it bounds a form of relative error in 

vector. This is crucial in the rare event context where the 

Notice for example that removing the ij (which are all non-n

inequality shows that the relative error in the invariant prob

errors in the non-zero entries of FIt is not dicult to constr

If and with for some

X

j6=i

Fij 6=0

Pi[tj < ti]

1

⇠

L

# of MCMC steps per window

i

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k)

zT= zTF

z ̃T =  ̃zTF ̃

Fij , F ̃ij  Sij

1

2(f) ⇠

L

# of MCMC steps per window

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k)

zT= zTF

z ̃T =  ̃zTF ̃

Fij , F ̃ij  Sij

1

Ri

jk =

i

jk

pvar⇡i ( j )

pvar⇡i ( k)

zT= zTF

z ̃T =  ̃zTF ̃

Fij , F ̃ij  Sij

1

sub-stochastic matrix then S

Fik

tr(Ri) ⇠

1

# of MCMC steps

Pi[tj < ti]  Fij

Fi,i±1 =

1

2

1

1 + exp ⇥LV  i±1

L



 V  i

L

⇤ 

1

2

1

1 + e

where

and and are respectively the first time to arrive at state and absorption 

time for the absorbing Markov process governed by 

⌧i

⌧

⌧i

⌧

Fik

e

F

= X

i=1

zi he

W

ii

eF= E

h

eW(t)

i

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

i

j

S

⌧

A. Dinner, E. Thiede, B. Van Koten, and JW, SIMAX 36(3), 917-941 (2015)
A typical US application: A free energy surface for the 

 alanine dipeptide 

184114-6 Tempkin et al. J. Chem. Phys. 140, 184114 (2014)

0 10 20 30 40 50 60

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

iteration

error

FG only

Multilevel with bad CG model

Multilevel with good CG model

Bad CG only

Good CG only

FIG. 6. String method for the Müller-Brown model. We plot the error (as

computed by Eq. (14)) at each iteration of the indicated string methods: the

string method with only the FG model (◦), the multilevel string method with

bad CG model (×), the multilevel string method with good CG model (+),

the string method applied directly to the bad CG model (∗), and the string

method applied directly to the good CG model (!).

smoothed FG model (i.e., the fixed point of S), the agreement

is not perfect. The two paths are plotted over isoenergetic contours of the good CG potential in the top panel of Fig. 5. The

difference (as measured by Eq. (14)) between the fixed point

N

H

H

N

H3C CH3

O

O

H3C H

FIG. 7. The collective variables used in the string method are the φ and ψ

backbone dihedral angles (above). The structures used for the PNM basins

are the c7ax and c7eq configurations (below).

used to constrain bonds involving hydrogen atoms.34 The

FG calculations are performed in NAMD2.9,35 and the CG

calculations are performed using the PNM module in

CHARMM (version c35b3).1

The collective variables (CVs) for the string procedure

are the φ and ψ backbone dihedrals (Figure 7). The string has

16 images. The initial position of the string is obtained by linearly interpolating between the points (φ, ψ) = (70.5, –69.5)

and (φψ) = (–827735)We generate initial structures
emcee

US (temp)

US (cv)

1.0

0.8

0.6

0.4

0.2

0

1.0

0.8

0.6

0.4

0.2

0

1.0

0.8

0.6

0.4

0.2

0

0 0.2 0.4 0.6 0.8 1.0

Joint posterior of cosmological values 

J = j

⌦m

Z

f(x)⇡(dx) = X

i

zi hfii

"

t=0

#

i

J = j

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

: mean dimensionless matter

: vacuum energy density

US with two different choices of bias functions:

(tempering) 

hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

⇡i / (⇡)

↵i

↵i 2 (0, 1]

⌦

Z

f(x)⇡(dx) = X

i

zi hfii

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi h

⇡i / (⇡)

↵i

↵i 2 (0, 1]

⌦

Z

f(x)⇡(dx) = X

i

zi hfi

(collective variable) bias functions restrict along strips 

parallel with long axis of level sets

In both cases US gives much more accurate representation 

 of the tails of the posterior 

In particular we obtain a more accurate estimate of the 

probability that the expansion of the universe is decelerating, 

contradicting previous (higher) estimates.

C. Matthews, JW, A. Kravstov, and E. Jennings, in preparation
new applications + improvements for current applications

E. Thiede, B. Van Koten, A. Dinner, and JW, J. Chem. Phys. 145(8), 084115 (2016)
We estimate expectations of the 

form:

Associate with an index 

process 

that labels regions in both 

space and time:

t t = ⌧

X(t)

J = j

t = 0

Trajectory Stratification Framework

zT= zTF, X

i

zi = 1

E

"⌧

X1

t=0

f(t, X(t))

#

X

i

 i = 1

 log i

Z

f(x)⇡(dx) = X

i

zi hfii

dXt = rV (Xt) dt + p2" dWt

V

J(t) 2 {0, 1, 2,...,n}

"1

{zi}

zT= zTF, X

i

zi = 1

E

"⌧

X1

t=0

f(t, X(t))

#

X

i

 i = 1

dXt = rV (Xt) dt + p2" dW

X(t)

J(t) 2 {0, 1, 2,...,n}

"1

{zi}

zT= zTF, X

i

zi = 1

E

"⌧

X1

t=0

f(t, X(t))

#

X

 i = 1

2 {0, 1, 2,...,n}

"1

⌧

zTF, X

i

zi = 1

(t, X(t))

#

= X

i

zi hfii

X

i

 i = 1

is a first hitting time for 

X(t)

J(t) 2 {0, 1, 2,...,n}

"1

(t, X(t))

zT= zTF, X

i

zi = 1

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

X

i

 i = 1

 log i

Z

f(x)⇡(dx) = Xzhfi

J. O. B. Tempkin, B. Van Koten, J. C. Mattingly, A. R. Dinner, J. Weare. arXiv:1610.09426
 

kcal/mol

Set A Set B

184114-6 Tempkin et al. J. Chem. Phys. 140, 184114 (2014)

0 10 20 30 40 50 60

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

iteration

error

FG only

Multilevel with bad CG model

Multilevel with good CG model

Bad CG only

Good CG only

FIG. 6. String method for the Müller-Brown model. We plot the error (as

computed by Eq. (14)) at each iteration of the indicated string methods: the

string method with only the FG model (◦), the multilevel string method with

bad CG model (×), the multilevel string method with good CG model (+),

the string method applied directly to the bad CG model (∗), and the string

method applied directly to the good CG model (!).

N

H

H

N

H3C CH3

O

O

H3C H

FIG. 7. The collective variables used in the string method are the φ and ψ

backbone dihedral angles (above). The structures used for the PNM basins

are the c7ax and c7eq configurations (below).

used to constrain bonds involving hydrogen atoms.34 The

FG calculations are performed in NAMD2.9,35 and the CG

calculations are performed using the PNM module in

Illustrative example: Computing set hitting probabilities in the 

isomerization of the alanine dipeptide

Compute the probability 

that the system hits the set B before hitting the set A and 

before time

"1

PBA = Px[⌧B < min{⌧A, ⌧max}]

ai = P[J(0) = i]

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

 log i

J(t) 2 {0, 1, 2,...,n}

"1

⌧max

ai = P[J(0) = i]

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j
A

B

PBA ⇡

NB

N

This probability can be very difficult to estimate when the 

probability of hitting B is very small.



X(0)

Illustrative example: Computing set hitting probabilities in the 

isomerization of the alanine dipeptide

(0) = 58.0

f(t, X(t)) = 1B(X(t))

g ypp

ility, PBA(⌧max), of hitting a set B before a separate set A and before a fixed time

max > 0 given that the system is at a point X(0) 

2/ A [ B at time t = 0. In the

se where X(0) and B are separated by a large free energy barrier while X(0) and A

e not, computing PBA(⌧max) can be challenging since trajectories that contribute

PBA(⌧max) are rare in direct simulations. To compute PBA(⌧max) via the scheme

subsection 3.3, we let the stopping time ⌧ be the minimum of ⌧max and the first

me, t, at which X(t1) is in either A or B, i.e., ⌧  1 = min{⌧A, ⌧B, ⌧max  1}

here ⌧A and ⌧B are the first times that X(t) enters the sets A and B respectively.

rictly speaking, to write ⌧ in the form in (9), we need to replace (t, X(t) ) in that

uation by (t, X(t1), X(t) ). The set D corresponding to our choice of ⌧ is then

= {(t, x, y) : t < ⌧max,x /2 (A [ B)}. As we have already mentioned, this can be

ne without further modification of the scheme. Then f(t, X(t) ) in (10) is

6) f(t, X(t)) = 1B(X(t) ).

The system that we simulate is the alanine dipeptide (CH3-CONH-C↵H(CH3)-

This manuscript is for review purposes only.
To compute using the non-stationary trajectory stratification, we 

divide the trajectory space along the phi dihedral and along time:

f(t, X(t)) = 1B(X(t))

E

"⌧

X1

t=0

f(t, X(t))#= X

n

i=1

zihfii

A

B



Illustrative example: Computing set hitting probabilities in the 

isomerization of the alanine dipeptide

"1

PBA

ai = P[J(0) = i]

zT= zTG + aT

E

"⌧

X1

t=0

f(t, X(t))

#

= X

i

zi hfii

J = j

 log i

m to 1. We choose t

start and t

end 

to equally divide the interval [0, ⌧ ), where, in our

xamples, ⌧ is a fixed time horizon. The function j is largest when t 2 [t

j

start, tj

end 

)

nd (x) = j . The supports of the various j correspond to products of overlapping

tervals in the  variable, but non-overlapping intervals in time. The fact that j

epends on time is essential in our examples.

5.2. Finite-Time Hitting Probability. In this section we compute the prob

bility, PBA(⌧max), of hitting a set B before a separate set A and before a fixed time

max > 0 given that the system is at a point X(0) 

2/ A [ B at time t = 0. In the

ase where X(0) and B are separated by a large free energy barrier while X(0) and A

e not, computing PBA(⌧max) can be challenging since trajectories that contribute

PBA(⌧max) are rare in direct simulations. To compute PBA(⌧max) via the scheme

subsection 3.3, we let the stopping time ⌧ be the minimum of ⌧max and the first

me, t, at which X(t1) is in either A or B, i.e., ⌧  1 = min{⌧A, ⌧B, ⌧max  1}

here ⌧A and ⌧B are the first times that X(t) enters the sets A and B respectively.

rictly speaking, to write ⌧ in the form in (9), we need to replace (t, X(t) ) in that

quation by (t, X(t1), X(t) ). The set D corresponding to our choice of ⌧ is then

= {(t, x, y) : t < ⌧max,x /2 (A [ B)}. As we have already mentioned, this can be

one without further modification of the scheme. Then f(t, X(t) ) in (10) is

6) f(t, X(t)) = 1B(X(t) ).

The system that we simulate is the alanine dipeptide (CH3-CONH-C↵H(C 

H3)-
Repeated direct integration (106 independent simulations) 

NEUS stratification

zi

zi

— Direct 

— Stratification
Direct Integration NEUS Stratification zi zi

Choosing close to A makes estimating 

 a hard rare event problem

(0) = 91.0

PBA

i = P[J(0) = i]

zT= zTG + aT

(t, X(t))

#

= X

i

zi hfii

J = j

 log i
Software Tools: Enhanced Sampling Toolkit

I’ve developed a set of software tools designed to facilitate rapid 

prototyping of enhanced sampling algorithms. 

MD Engines: 

• LAMMPS, OpenMM, CHARMM, 

GROMACS (potential future 

implementation)

• Dynamics are executed in 

widely available and popular MD 

codes

• Computationally expensive part 

of the code is still executed in 

fast MD codes

Enhanced Sampling Toolkit

Enhanced Sampling Algorithms: 

• Allows developers to write 

algorithms in 100% Python 

• Enables rapid prototyping in 

developer friendly languages 

• Enables clean and concise 

writing of algorithmic code