ADAPTIVE TWO-STAGE INTEGRATORS FOR SAMPLING ALGORITHMS BASED ON - - PowerPoint PPT Presentation

• E. Akhmatskayaa,c, M. Fernández-Pendása, T. Radivojevića, J. M. Sanz-Sernab

aBasque Center for Applied Mathematics (BCAM), Bilbao, Spain bDepartamento de Matemáticas, Universidad Carlos III de Madrid, Madrid, Spain cIKERBASQUE, Basque Foundation for Science, Bilbao, Spain

ADAPTIVE TWO-STAGE INTEGRATORS FOR SAMPLING ALGORITHMS BASED ON HAMILTONIAN DYNAMICS

Let is Hamiltonian Equations of motion associated with H:

θ – position, p – momentum, U – potential energy, M – mass matrix, , D - dimension (θ,p)∈R2D

SAMPLING WITH HAMILTONIAN DYNAMICS

H(θ,p)= 1

2pTM−1p+U(θ)

The algorithms sampling with Hamiltonian dynamics: ª Can be viewed as special cases of Generalized Hybrid Monte Carlo (GHMC) [Kennedy, Pendleton, 2001] ª Numerically solve Hamiltonian equations

dθ dt =M−1p; dp dt =-Uθ(θ)

GHMC AT GLANCE

Initialization: positions θ, momenta p (drawn from Gaussian distribution), step size h, length of trajectories L, number of samples N, φ Proposal: numerical integration of Hamiltonian equations with the symplectic method Ψh over L steps and step-size h: Sampling: Metropolis accept / reject test Momentum refreshment:

ΨT (θ,p) = ( ! θ , ! p ),T = Lh

( !

θ , ! p ) = FΨT (θ,p) with probability min(1,exp(-βΔH)) (θ,p)

• therwise

# $ % % % & % % % ΔH = H(ΨT (θ,p))− H(θ,p) = H( FΨT (θ,p))− H(θ,p)

Momentum flip F : (θ,p) = (θ,-p) ! p ! u " # $ % & ' = cosϕ sinϕ −sinϕ cosϕ " # $ $ % & ' ' F p u " # $ $ % & ' ', 0 <ϕ ≤ π / 2, u ~ N(0,β −1M)

METHODS & APPLICATIONS

Molecular Simulations Computational Statistics

MC

GHMC GSHMC MMHMC

Sample with Shadow Hamiltonians

HMC HMC

Hybrid Hamiltonian φ=π/2

METHODS & APPLICATIONS

Molecular Simulations Computational Statistics

MC

GHMC GSHMC MMHMC

Sample with Shadow Hamiltonians

HMC HMC

Hybrid Hamiltonian φ=π/2 Adaptive parameters

NUTS

Shadow Hamiltonians

SHMC / S2HMC

Delayed rejections

LAHMC / ECGHMC

IN FOCUS

Molecular simulations: ª Hybrid Monte Carlo (HMC)[Duane, et. al, 1987] ª Generalized shadow hybrid Monte Carlo (GSHMC)

[Akhmatskaya, Reich, 2008]

HMC: GHMC with a complete momentum update GSHMC: Importance sampling GHMC

IN FOCUS

Computational statistics: ª Hamiltonian Monte Carlo (HMC)[Neil, 1993] ª Mixed & Match Hamiltonian Monte Carlo (MMHMC)

[Radivojevic, Akhmatskaya, 2016]

HMC: Hybrid Monte Carlo formulated for statistics MMHMC: GSHMC adjusted for statistics ª For the target density π(θ) of position vector θ, momenta

p conjugate to θ, with a ‘mass’ matrix M (a preconditioner)

Potential function: Hamiltonian:

U(θ) = −log(π(θ)) H(θ,p)= 1

2pTM−1p+U(θ)

WHAT INTEGRATOR TO CHOOSE?

Other suggestions?

NUMERICAL INTEGRATORS FOR SIMULATING HAMILTONIAN DYNAMICS

Wish list: ª Reversible ª Symplectic (preserve the symplectic structure of the Hamiltonian dynamics) => volume preserving ª Low integration error ª Computationally efficient ª Stable

VERLET INTEGRATOR

Verlet or Verlet / Stormer or leapfrog: a golden standard for Hamiltonian dynamics based simulation methods ª Reversible & symplectic ª Second order of accuracy (global error of O (h2), h – time step) ª One force evaluation per step ª Stability interval: (0,2) ª Easy to implement ª Verlet as a splitting integrator:

• t-flow of X={A, B}

Can we beat it?

ψh =ϕh/2

B ϕh A ϕh/2 B

ϕt

X

SPLITTING INTEGRATORS

Potential competitors to Verlet are members of one-parameter family of 2-stage splitting integrators of Hamiltonian system with

b – is a parameter of the family; b=0.25 for Verlet

2-stage splitting integrators: ª Capable to provide higher accuracy than Verlet (substeps h/2) ª Can achieve bigger time steps (stability interval up to: (0,4)) but ª Two force evaluations per step ª Fair comparison with Verlet: equal numbers of force evaluations

ψh =ϕbh

B ϕh/2 A ϕ(1−2b)h B

ϕh/2

A ϕbh B

hVerlet = h2−stage / 2; LVerlet = 2L2−stage H(θ,p)= 1

2pTM−1p+U(θ)≡ A+B

2-STAGE INTEGRATORS: SPECIAL CASE 1

Minimum-error (ME) integrator by McLachlan (1995): ª Criterion: to minimize the error εper step ª Resulting parameter: b ≈ 0.1932 ª Performance:

ü outperforms Verlet at small time steps ü performance degrades for large h: stability interval (0, 2.55)

ε ≈ Ch

3 + O(h 5 );

h → 0

2-STAGE INTEGRATORS: SPECIAL CASE 2

BCSS integrator by Blanes, Casas and Sanz-Serna (2014): ª Criterion: to minimize expectation of the energy error, , for ª Resulting parameter: b ≈ 0.2113 ª Performance: ü Shows the best performance around h=2 ü Performance may drop for smaller or larger step sizes Any better values of b?

Δ = H(ψh,L(θ,p))− H(θ,p) 0 ≤ E(Δ)≤ ρ(h,b) ρ(h,b) = h4(2b2(1/ 2−b)h2 + 4b2 −6b+1)2 8(2−bh2)(2−(1/ 2−b)h2)(1−b(1/ 2−b)h2) E(Δ) 0 < h < 2

ADAPTIVE INTEGRATION APPROACH AIA

Idea: to present a parameter b as a function of internal properties (frequencies) of a simulated system and a simulation time step, i.e. , in such a way that: Given a simulation problem (ω)and a simulation time (Δt), b defines the unique integrator which guarantees the best conservation of energy for harmonic forces. ª The resulting integrator: a problem specific two-stage integrator ª Criterion: to minimize ª AIA uses the upper bound suggested in BCSS method

b = b(Δt,ω) max

0<h<h ρ(h,b), h ≡ h(Δt,ω)

Given a simulation problem (ω)and a simulation time (Δt): 1. Compute the non-dimensional quantity

• highest frequencies and dimension of the system,
• average energy error obtained from warm-up simulation with

2. If h>4 abort the integration 3. Find the optimal value of the parameter b as Molecular Simulation: S = 2

AIA: ALGORITHM

b = arg min

b∈(0,0.25) max h∈(0,h )ρ(h,b)

h = S  ωΔt  ω,D dH

VV

ΔtVV

(Schlick, 2002; Skeel, 1999)

Statistics: S = 2 dH

VV / 2

ωi

6 i=1 D

∑

" # $ $ % & ' '

6

ΔtVV

AIA: REMARKS

ª No computational overheads in simulations ª Extended to constrained dynamics: ü Presented the two-stage integrator as two concatenated velocity Verlet steps and combined with RATTLE

AIA FOR MODIFIED HAMILTONIAN METHODS

ª MAIA: AIA for methods sampling with modified Hamiltonian ª The same ideas as in AIA but the expected value of modified energy w.r.t. to the modified density is minimised ª MAIA requires: ü Modified Hamiltonians for 2-stage splitting integrators ü Upper bound for expectation of modified energies

Δ =  H(ψh,L(θ,p))−  H(θ,p)

• numerical time derivative of
• parameter of a 2-stage integrator

MODIFIED HAMILTONIANS FOR SPLITTING INTEGRATORS

4-th order modified Hamiltonians for 2-stage splitting integrators were derived in terms of quantities available during simulation

Model problem: Standard harmonic oscillator in non-dimensional variables (standard univariate Gaussian). The solution flow is given by depends on a integrator. Defining obtain:

MAIA: UNDERLYING ANALYSIS

ü Numerical solution is exact solution of modified Hamiltonian ü Numerical solution stays on ellipse ü Maximum Δ: numerical solution starts at ellipse vertices on major semi- axis

UPPER BOUND FOR EXPECTATION OF MODIFIED ENERGY ERROR

D-variate Gaussian target distribution (d coupled linear oscillators): ωj - the angular frequencies of the oscillators

E(Δ) ≤ ρ(ω jh)

j=1 D

∑

1 2 3 4 h 0.17 0.19 0.21 0.23 0.25 0.27 b

INTEGRATOR’S PARAMETER VS. TIME STEP

1 2 3 4 h 0.17 0.19 0.21 0.23 0.25 0.27 b

VV AIA BCSS ME

VV MAIA M-BCSS M-ME

Sampling with Hamiltonians Sampling with shadow Hamiltonians

EXPECTED ERROR VS. TIME STEP

1 2 3 4 h 0.1 0.2 0.3 0.4 Expected error 1 2 3 4 h 0.1 0.2 0.3 0.4 Expected error

VV AIA BCSS ME

VV MAIA M-BCSS M-ME

Sampling with Hamiltonians Sampling with shadow Hamiltonians

EXPECTED ERROR VS. TIME STEP: ZOOM

1 2 3 4 h 0.01 0.02 0.03 0.04 0.05 Expected error

VV AIA BCSS ME

VV MAIA M-BCSS M-ME

TIME STEP vs. EXPECTED ERROR

0.1 0.2 0.3 0.4 Expected error 1 2 3 4 h 0.1 0.2 0.3 0.4 Expected error 1 2 3 4 h

VV AIA BCSS ME

Sampling with Hamiltonians Sampling with shadow Hamiltonians

VV MAIA M-BCSS M-ME

0.01 0.02 0.03 0.04 0.05 Expected error 1 2 3 4 h

TIME STEP vs. EXPECTED ERROR: ZOOM

VV AIA BCSS ME

VV MAIA M-BCSS M-ME

(M)AIA IMPLEMENTATION

ª Can be implemented in any MD/HMC software with no

• verheads

ª Two major implementation steps: 1. Introduce AIA in a pre-processing unit – to be run once before a simulation

ü find b and pass to a simulator

2. Present 2-stage scheme in a kick/drift factorization form in a simulation module ª Current implementation ü MS: in multiHMC-GROMACS – BCAM in-house modified version of GROMACS 4.5.4 ü CS: in BCAM software package Hamiltonians in Computational Statistics (HaiCS)

TESTING (M)AIA

Molecular Simulation Spider toxin in membrane/water environment: Coarse grained unconstrained system with 7810 particles Villiin: Atomistic system with 9389 atoms. Hydrogens are constrained using the SHAKE/RATTLE algorithm Computational Statistics Gaussian distribution: D=1000, 2000 Metrics ESS – time normalized effective sample size

• Rel. ESS – relative ESS w.r.t. ESS achieved with Verlet

IACF – integrated autocorrelation function

2 3 4 5 Δt (fs) 10 20 30 40 50 60 70 80 90 100 AR (%) 2 3 4 5 Δt (fs) 10 20 30 40 50 60 70 80 90 100 AR (%) 20 25 30 35 40 45 50 Δt (fs) 10 20 30 40 50 60 70 80 90 100 AR (%)

ACCEPTANCE RATES: TOXIN & VILLIN

20 25 30 35 40 45 50 Δt (fs) 10 20 30 40 50 60 70 80 90 100 AR (%)

HMC GSHMC

VV AIA BCSS

VV MAIA M-BCSS

Toxin Villin

20 30 40 Δt (fs) 10 20 30 40 50 60 70 IACF x time (h)

IACF: TOXIN

20 30 40 Δt (fs) 10 20 30 40 50 60 70 IACF x time (h)

HMC

VV AIA BCSS

VV MAIA M-BCSS

GSHMC

TOXIN-BILAYER DISTANCE

1 2 3 4

Distance (nm)

0.5 1 1.5 2 2.5 3

True distribution VV BCSS AIA

5000 10000 15000 20000

Time (ps)

0.5 1 1.5 2 2.5 3

Distance (nm) VV BCSS AIA Target distance

HMC HMC

∆t

0.008 0.012 0.016 0.02 0.024

AR (%)

25 50 75 100

HMC ∆t

0.016 0.02 0.024 0.028

MMHMC ∆t

0.015 0.02 0.025 0.03 0.035

AR (%)

25 50 75 100

HMC ∆t

0.02 0.025 0.03 0.035 0.04

MMHMC

ACCEPTANCE RATES: GAUSSIAN

VV AIA BCSS ME VV MAIA M-BCSS M-ME

D=1000 D=2000

RELATIVE ESS: GAUSSIAN, D=1000

∆t

0.015 0.02 0.025 0.03 0.035

• rel. medESS

1 2 3 4 5 6 7

• rel. min ESS

1 2 3 4 5 6 7

HMC ∆t

0.02 0.025 0.03 0.035 0.04

MMHMC

VV AIA BCSS ME VV MAIA M-BCSS M-ME

RELATIVE ESS: GAUSSIAN, D=2000

∆t

0.008 0.012 0.016 0.02 0.024

• rel. medESS

2 4 6 8 10 12 14

• rel. min ESS

2 4 6 8 10 12 14

HMC ∆t

0.016 0.02 0.024 0.028

MMHMC

VV AIA BCSS ME VV MAIA M-BCSS M-ME

ESS: GAUSSIAN

h

1 1.5 2

min ESS

500 1000 1500

HMC h

1 1.5 2 2.5 500 1000 1500

MMHMC

VV AIA BCSS ME VV MAIA M-BCSS M-ME

h

1 1.5 2

min ESS

500 1000 1500 2000

HMC h

1.5 2 2.5 500 1000 1500 2000

MMHMC

D=1000 D=2000

CONCLUSIONS

ª We present an alternative to the standard Verlet integrator, (M)AIA, which offers, for any chosen simulation problem and step size, a system-specific two-stage splitting integrator to provide the best conservation of energy for harmonic forces ª The method was derived for sampling with Hamiltonians and modified Hamiltonians and extended to constrained dynamics ª Testing M(AIA) in molecular simulation and statistical applications demonstrated its superiority over velocity Verlet, BCSS and ME integrators ü never worse than VV, BCSS and ME by design ü leads to higher acceptance rates and better sampling ü allows for longer time steps ü no computational overheads in simulations

WORK IN PROGRESS

ª M(AIA) for n-stage splitting integrators (n>2) ª M(AIA) for multiple-time stepping methods ª Predicting a range of optimal time steps using the (M)AIA underlying analysis

Guggenheim, Bilbao BCAM

Modeling

• deling & Simulat

imulation ion in in Lif Life e & Mat ater erials ials Sciences ciences: : MSLM LMS