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

ADAPTIVE TWO-STAGE INTEGRATORS FOR SAMPLING ALGORITHMS BASED ON HAMILTONIAN DYNAMICS E. Akhmatskaya a,c , M. Fernández-Pendás a , T. Radivojevi ć a , J. M. Sanz-Serna b a Basque Center for Applied Mathematics (BCAM), Bilbao, Spain b Departamento de Matemáticas, Universidad Carlos III de Madrid, Madrid, Spain c IKERBASQUE, Basque Foundation for Science, Bilbao, Spain

SAMPLING WITH HAMILTONIAN DYNAMICS H ( θ , p) = 1 2 p T M − 1 p +U ( θ ) Let is Hamiltonian θ – position, p – momentum, U – potential energy, ( θ , p) ∈ R 2D M – mass matrix , , D - dimension Equations of motion associated with H : dp d θ d t =M − 1 p; d t = - 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

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 : p ), T = L h Ψ T ( θ , p ) = ( ! θ , ! Sampling: Metropolis accept / reject test # F Ψ T ( θ ,p ) with probability min(1,exp(- β Δ H )) % % % ( ! θ , ! p ) = $ % otherwise ( θ ,p ) % % & flip F : ( θ , p ) = ( θ ,- p ) Momentum Δ H = H ( Ψ T ( θ ,p )) − H ( θ ,p ) = H ( F Ψ T ( θ ,p )) − H ( θ ,p ) Momentum refreshment: " % " % " p % cos ϕ sin ϕ p ! ' F u ~ N (0, β − 1 M ) $ ' ' , 0 < ϕ ≤ π / 2, $ ' ' = $ $ $ u − sin ϕ cos ϕ u ! # & # & # &

METHODS & APPLICATIONS Molecular Simulations Computational Statistics GHMC HMC HMC GSHMC MMHMC Hybrid Hamiltonian Sample with Shadow φ = π /2 Hamiltonians MC

METHODS & APPLICATIONS Molecular Simulations Computational Statistics GHMC HMC HMC GSHMC MMHMC Hybrid Hamiltonian Sample with Shadow φ = π /2 Hamiltonians MC Delayed rejections Shadow Hamiltonians Adaptive parameters LAHMC / ECGHMC SHMC / S2HMC NUTS

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) U ( θ ) = − log ( π ( θ )) Potential function: H ( θ , p) = 1 2 p T M − 1 p +U ( θ ) Hamiltonian:

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 ( h 2 ), h – time step) ª One force evaluation per step ª Stability interval: (0,2) ª Easy to implement B  ϕ h A  ϕ h /2 B ψ h = ϕ h /2 ª Verlet as a splitting integrator: - t - flow of X ={ A, B } X ϕ t Can we beat it?

SPLITTING INTEGRATORS Potential competitors to Verlet are members of one-parameter family of 2-stage splitting integrators of Hamiltonian system with H ( θ , p) = 1 2 p T M − 1 p +U ( θ ) ≡ A + B B  ϕ h /2 A  ϕ (1 − 2 b ) h A  ϕ bh B B  ϕ h /2 ψ h = ϕ bh 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 Verlet = h 2 − stage / 2; L Verlet = 2L 2 − stage

2-STAGE INTEGRATORS: SPECIAL CASE 1 Minimum-error (ME) integrator by McLachlan (1995): ª Criterion: to minimize the error ε per step 5 ); 3 + O ( h ε ≈ Ch h → 0 ª Resulting parameter: b ≈ 0.1932 ª Performance: ü outperforms Verlet at small time steps ü performance degrades for large h : stability interval (0, 2.55)

2-STAGE INTEGRATORS: SPECIAL CASE 2 BCSS integrator by Blanes, Casas and Sanz-Serna (2014): ª Criterion: to minimize expectation of the energy error, , E( Δ ) for 0 < h < 2 Δ = H ( ψ h , L ( θ , p )) − H ( θ , p ) 0 ≤ E( Δ ) ≤ ρ ( h , b ) h 4 (2 b 2 (1/ 2 − b ) h 2 + 4 b 2 − 6 b + 1) 2 ρ ( h , b ) = 8(2 − bh 2 )(2 − (1/ 2 − b ) h 2 )(1 − b (1/ 2 − b ) h 2 ) ª 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 ?

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: b = b ( Δ t , ω ) 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 max 0 < h < h ρ ( h , b ), h ≡ h ( Δ t , ω ) ª Criterion: to minimize ª AIA uses the upper bound suggested in BCSS method

AIA: ALGORITHM Given a simulation problem ( ω )and a simulation time ( Δ t ): h = S  ω Δ t 1. Compute the non-dimensional quantity Molecular Simulation: S = 2 (Schlick, 2002; Skeel, 1999) " % D 6 d H VV / 2 ∑ $ ω i ' 6 $ ' # & i = 1 Statistics : S = 2 Δ t VV d H  ω , D - highest frequencies and dimension of the system, VV Δ t VV - average energy error obtained from warm-up simulation with 2. If h>4 abort the integration Find the optimal value of the parameter b as 3. b = arg min b ∈ (0,0.25) max h ∈ (0, h ) ρ ( h , b )

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 Δ =  H ( ψ h , L ( θ , p )) −  energy H ( θ , p ) 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

MODIFIED HAMILTONIANS FOR SPLITTING INTEGRATORS 4-th order modified Hamiltonians for 2-stage splitting integrators were derived in terms of quantities available during simulation - numerical time derivative of - parameter of a 2-stage integrator

MAIA: UNDERLYING ANALYSIS Model problem: Standard harmonic oscillator in non-dimensional variables (standard univariate Gaussian). The solution flow is given by depends on a integrator. Defining obtain: ü 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): D E ( Δ ) ≤ ∑ ρ ( ω j h ) j = 1 ω j - the angular frequencies of the oscillators

INTEGRATOR’S PARAMETER VS. TIME STEP Sampling with Hamiltonians Sampling with shadow Hamiltonians 0.27 0.27 0.25 0.25 0.23 0.23 b b 0.21 0.21 0.19 0.19 0.17 0.17 0 1 2 3 4 0 1 2 3 4 h h VV VV AIA MAIA BCSS M-BCSS ME M-ME

EXPECTED ERROR VS. TIME STEP Sampling with Hamiltonians Sampling with shadow Hamiltonians 0.4 0.4 VV VV AIA MAIA M-BCSS BCSS 0.3 0.3 Expected error Expected error M-ME ME 0.2 0.2 0.1 0.1 0 0 0 1 2 3 4 0 1 2 3 4 h h

EXPECTED ERROR VS. TIME STEP: ZOOM 0.05 VV AIA BCSS 0.04 ME Expected error 0.03 0.02 VV 0.01 MAIA M-BCSS M-ME 0 0 1 2 3 4 h

TIME STEP vs. EXPECTED ERROR Sampling with Hamiltonians Sampling with shadow Hamiltonians 4 4 3 3 2 h 2 h VV VV 1 AIA 1 MAIA BCSS M-BCSS ME M-ME 0 0 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 Expected error Expected error

TIME STEP vs. EXPECTED ERROR: ZOOM 4 3 2 h 1 VV VV MAIA AIA M-BCSS BCSS M-ME ME 0 0 0.01 0.02 0.03 0.04 0.05 Expected error