Computing ergodic limits for SDEs M.V. Tretyakov School of - PowerPoint PPT Presentation

Computing ergodic limits for SDEs M.V. Tretyakov School of Mathematical Sciences, University of Nottingham, UK Talk at the workshop “Stochastic numerical algorithms, multiscale modelling and high-dimensional data analytics”, ICERM, 21st July 2016

Plan of the talk t � 1 � ϕ erg = ϕ ( x ) ρ ( x ) dx = lim t →∞ E ϕ ( X ( t )) = lim ϕ ( X ( s )) ds a . s . t t →∞ 0

Plan of the talk t � 1 � ϕ erg = ϕ ( x ) ρ ( x ) dx = lim t →∞ E ϕ ( X ( t )) = lim ϕ ( X ( s )) ds a . s . t t →∞ 0 Introduction Examples of Langevin-type equations and stochastic gradient systems Geometric integrators for Langevin equations and the gradient system ‘Non-Markovian’ scheme for stochastic gradient systems [Davidchack, Ouldridge&T. J Chem Phys 2015] [Leimkuhler, Matthews,T 2014]

Introduction Hamiltonian H ( r , p ) canonical ensemble ( NVT ) ρ ( r , p ) ∝ exp( − β H ( r , p )) , where β = 1 / ( k B T ) > 0 is an inverse temperature. t 1 � � ϕ erg = ϕ ( r ) ρ ( r , p ) drdp = lim t →∞ E ϕ ( R ( t )) = lim ϕ ( R ( s )) ds a . s . t t →∞ 0

Rigid Body Dynamics Following [Miller III et al J. Chem. Phys. , 2002] n 3 H ( r , p , q , π ) = p T p � � V l ( q j , π j ) + U ( r , q ) , 2 m + (1) j =1 l =1 T T T T ) T ∈ R 3 n , and p = ( p 1 ) T ∈ R 3 n are the where r = ( r 1 , . . . , r n , . . . , p n T T ) T ∈ R 4 n , center-of-mass coordinates and momenta; q = ( q 1 , . . . , q n q j = ( q j 3 ) T ∈ R 4 , | q j | = 1 , are the rotational coordinates in the 0 , q j 1 , q j 2 , q j T T ) T ∈ R 4 n , are the angular quaternion representation, π = ( π 1 , . . . , π n momenta; 3 3 V l ( q , π ) = 1 1 � 2 = 1 � � π T S l q 8 π T S ( q ) DS T ( q ) π, q , π ∈ R 4 , l = 1 , 2 , 3 , � 8 I l l =1 l =1 I l – the principal moments of inertia and the constant 4-by-4 matrices S l : ( − q 1 , q 0 , q 3 , − q 2 ) T , S 2 q = ( − q 2 , − q 3 , q 0 , q 1 ) T , S 1 q = ( − q 3 , q 2 , − q 1 , q 0 ) T , S 0 = diag(1 , 1 , 1 , 1) T . S 3 q = S ( q ) = [ S 0 q , S 1 q , S 2 q , S 3 q ] , D = diag(0 , 1 / I 1 , 1 / I 2 , 1 / I 3 ) .

Langevin thermostat for Rigid Body Dynamics P j dR j R j (0) = r j , = m dt , (2) � 2 m γ dP j f j ( R , Q ) dt − γ P j dt + dw j ( t ) , P j (0) = p j , = β 1 dQ j 4 S ( Q j ) DS T ( Q j )Π j dt , Q j (0) = q j , | q j | = 1 , = (3) 3 1 1 � d Π j Π j T S l Q j � S l Π j dt + F j ( R , Q ) dt − Γ J ( Q j )Π j dt � = 4 I l l =1 � 3 2 M Γ Π j (0) = π j , q j T π j = 0 , j = 1 , . . . , n , � S l Q j dW j + l ( t ) , β l =1 where f j ( r , q ) = −∇ r j U ( r , q ) ∈ R 3 , F j ( r , q ) = − ˜ ∇ q j U ( r , q ) ∈ T q j S 3 , ( w T , W T ) T = ( w 1 T , . . . , w n T , W 1 T , . . . , W n T ) T is a (3 n + 3 n )-dimensional standard Wiener process; γ ≥ 0 and Γ ≥ 0 are the friction coefficients for the translational and rotational motions, and 3 � J ( q ) = MS ( q ) DS T ( q ) / 4 , M = 4 / 1 / I l . l =1 Davidchack, Ouldridge&T. J Chem Phys 2015

Langevin-type equations The Ito interpretation of the SDEs (2)–(3) coincides with its Stratonovich interpretation. The solution of (2)–(3) preserves the quaternion length | Q j ( t ) | = 1 , j = 1 , . . . , n , for all t ≥ 0 . (4) The solution of (2)–(3) automatically preserves the constraint: Q j T ( t )Π j ( t ) = 0 , j = 1 , . . . , n , for t ≥ 0 (5) Assume that the solution X ( t ) = ( R T ( t ) , P T ( t ) , Q T ( t ) , Π T ( t )) T of (2)–(3) is an ergodic process on { x = ( r T , p T , q T , π T ) T ∈ R 14 n : D = q j T π j = 0 , | q j | = 1 , j = 1 , . . . , n } . Then it can be shown that the invariant measure of X ( t ) is Gibbsian with the density ρ ( r , p , q , π ) on D : ρ ( r , p , q , π ) ∝ exp( − β H ( r , p , q , π )) (6) Davidchack, Ouldridge&T. J Chem Phys 2015

Langevin equations and quasi-symplectic integrators P j dR j R j (0) = r j , = m dt , (9) � 2 m γ dP j f j ( R , Q ) dt − γ P j dt + dw j ( t ) , P j (0) = p j , = β 1 dQ j 4 S ( Q j ) DS T ( Q j )Π j dt , Q j (0) = q j , | q j | = 1 , = (10) 3 1 1 � d Π j � Π j T S l Q j � S l Π j dt + F j ( R , Q ) dt − Γ J ( Q j )Π j dt = 4 I l l =1 � 3 2 M Γ Π j (0) = π j , q j T π j = 0 , j = 1 , . . . , n , � S l Q j dW j + l ( t ) , β l =1

Langevin equations and quasi-symplectic integrators P j dR j R j (0) = r j , = m dt , (9) � 2 m γ dP j f j ( R , Q ) dt − γ P j dt + dw j ( t ) , P j (0) = p j , = β 1 dQ j 4 S ( Q j ) DS T ( Q j )Π j dt , Q j (0) = q j , | q j | = 1 , = (10) 3 1 1 � d Π j � Π j T S l Q j � S l Π j dt + F j ( R , Q ) dt − Γ J ( Q j )Π j dt = 4 I l l =1 � 3 2 M Γ Π j (0) = π j , q j T π j = 0 , j = 1 , . . . , n , � S l Q j dW j + l ( t ) , β l =1 Let D 0 ∈ R d , d = 14 n , be a domain with finite volume. The transformation x = ( r , p , q , π ) �→ X ( t ) = X ( t ; x ) = ( R ( t ; x ) , P ( t ; x ) , Q ( t ; x ) , Π( t ; x )) maps D 0 into the domain D t .

Langevin equations and quasi-symplectic integrators � dX 1 . . . dX d V t = (7) D t � D ( X 1 , . . . , X d ) � � � dx 1 . . . dx d . � � = � � D ( x 1 , . . . , x d ) � D 0 The Jacobian J is equal to J = D ( X 1 , . . . , X d ) D ( x 1 , . . . , x d ) = exp ( − n (3 γ + Γ) · t ) . (8) [Milstein&T. IMA J. Numer. Anal. 2003 (also Springer 2004)], [Davidchack, Ouldridge&T. J Chem Phys 2015]

The stochastic gradient system It is easy to verify that � exp( − β H ( r , p , q , π )) d p d π (9) D mom exp( − β U ( r , q )) =: ˜ ρ ( r , q ) , ∝ where ( r T , q T ) T ∈ D ′ = { ( r T , q T ) T ∈ R 7 n : | q j | = 1 } and the domain of T ) T ∈ R 7 n : q T π = 0 } . conjugate momenta D mom = { ( p T , π

The stochastic gradient system It is easy to verify that � exp( − β H ( r , p , q , π )) d p d π (9) D mom exp( − β U ( r , q )) =: ˜ ρ ( r , q ) , ∝ where ( r T , q T ) T ∈ D ′ = { ( r T , q T ) T ∈ R 7 n : | q j | = 1 } and the domain of T ) T ∈ R 7 n : q T π = 0 } . conjugate momenta D mom = { ( p T , π We introduce the gradient system in the form of Stratonovich SDEs: � d R = υ 2 υ m f ( R , Q ) dt + m β d w ( t ) , R (0) = r , (10) � 3 dQ j = Υ 2Υ S l Q j ⋆ dW j � M F j ( R , Q ) dt + l ( t ) , (11) M β l =1 Q j (0) = q j , | q j | = 1 , j = 1 , . . . , n , where “ ⋆ ” indicates the Stratonovich form of the SDEs, parameters υ > 0 and Υ > 0 control the speed of evolution of the gradient system (10)–(11), f = ( f 1 T , . . . , f n T ) T and the rest of the notation is as in (2)–(3). [Davidchack, Ouldridge&T. J Chem Phys 2015]

The gradient thermostat for rigid body dynamics This new gradient thermostat possesses the following properties. As in the case of (2)–(3), the solution of (10)–(11) preserves the quaternion length (4). Assume that the solution X ( t ) = ( R T ( t ) , Q T ( t )) T ∈ D ′ of (10)–(11) is an ergodic process. Then, by the usual means of the stationary Fokker-Planck equation, one can show that its invariant measure is Gibbsian with the density ˜ ρ ( r , q ) from (9).

Langevin integrators Davidchack, Ouldridge&T. J Chem Phys 2015 For simplicity we use a uniform time discretization of a time interval [0 , T ] with the step h = T / N . Goal: to construct integrators quasi-symplectic preserve | ¯ Q j ( t k ) | = 1 , j = 1 , . . . , n , for all t ≥ 0 automatically preserve ¯ Q j T ( t k )¯ Π j ( t k ) = 0 , j = 1 , . . . , n , for t ≥ 0 automatically of weak order 2 To this end: stochastic numerics+splitting techniques [see e.g. Milstein&T, Springer 2004] the deterministic symplectic integrator from [Miller III et al J. Chem. Phys. , 2002]

‘Langevin A’ integrator The first integrator is based on splitting the Langevin system in P j dR j R j (0) = r j , = m dt , (12) � 2 m γ dP j f j ( R , Q ) dt + dw j ( t ) , = β 1 dQ j 4 S ( Q j ) DS T ( Q j )Π j dt , = (13) 3 1 1 � Π j T S l Q j � d Π j � S l Π j dt + F j ( R , Q ) dt = 4 I l l =1 � 3 2 M Γ � S l Q j dW j + l ( t ) , j = 1 , . . . , n , β l =1

‘Langevin A’ integrator The first integrator is based on splitting the Langevin system in P j dR j R j (0) = r j , = m dt , (12) � 2 m γ dP j f j ( R , Q ) dt + dw j ( t ) , = β 1 dQ j 4 S ( Q j ) DS T ( Q j )Π j dt , = (13) 3 1 1 � Π j T S l Q j � d Π j � S l Π j dt + F j ( R , Q ) dt = 4 I l l =1 � 3 2 M Γ � S l Q j dW j + l ( t ) , j = 1 , . . . , n , β l =1 and the deterministic system of linear differential equations π j = − Γ J ( q j ) π j , j = 1 , . . . , n . p = − γ p , ˙ ˙ (14)

‘Langevin B’ integrator � 2 m γ d P I = − γ P I dt + d w ( t ) , β (15) � 3 2 M Γ � d Π j I = − Γ J ( q )Π j S l qdW j I dt + l ( t ); β l =1 P II II = 1 m dt , d P II = f ( R II , Q II ) dt , dQ j 4 S ( Q j II ) DS T ( Q j II )Π j d R II = II dt , (16) 3 F j ( R II , Q II ) dt + 1 1 � � d Π j � (Π j II ) T S l Q j S l Π j = II dt , j = 1 , . . . , n . II II 4 I l l =1