Monte Carlo in the Physical and Biological Sciences: Some Problems of Interest and Algorithms Paul Dupuis Division of Applied Mathematics Brown University ICERM 26 October 2012 Paul Dupuis (Brown University) 26 October 2012
Some problems of interest The ergodic problem—computing expected values with respect to a 1 stationary distribution (usually a Gibbs measure). Transition rate problems—computing the probability of transitions over 2 [ 0 , T ] , exit locations, and mean exit times with respect to metastable states [Jonathan Weare tutorial, Tuesday morning]. Paul Dupuis (Brown University) 26 October 2012
Examples A first example. Compute the average potential energy and other functionals with respect to a Gibbs measure of the form � π ( dx ) = e − V ( x ) /τ dx Z ( τ ) , and V is the potential of a (relatively) complex physical system. Paul Dupuis (Brown University) 26 October 2012
Examples A first example. Compute the average potential energy and other functionals with respect to a Gibbs measure of the form � π ( dx ) = e − V ( x ) /τ dx Z ( τ ) , and V is the potential of a (relatively) complex physical system. Primary interest is as the marginal (independent) distribution on spatial variables of stationary distribution of a Hamiltonian system. E.g., p 2 � n π ( dx , dp ) ∝ e − 1 τ V ( x ) − 1 j 2 m dxdp . ¯ j = 1 τ Paul Dupuis (Brown University) 26 October 2012
Examples A first example. Compute the average potential energy and other functionals with respect to a Gibbs measure of the form � π ( dx ) = e − V ( x ) /τ dx Z ( τ ) , and V is the potential of a (relatively) complex physical system. Primary interest is as the marginal (independent) distribution on spatial variables of stationary distribution of a Hamiltonian system. E.g., p 2 � n π ( dx , dp ) ∝ e − 1 τ V ( x ) − 1 j 2 m dxdp . ¯ j = 1 τ In this case the dynamical model is deterministic, and randomness only enters through initial condition x j ( t ) = 1 p j ( t ) = − ∂ ˙ mp j ( t ) , ˙ V ( x ( t )) . ∂ x j Paul Dupuis (Brown University) 26 October 2012
Examples One can use that π ( dx ) is the stationary distribution of the solution to √ dX = −∇ V ( X ) dt + 2 τ dW , as well as a variety of related discrete time models (will mention specific examples later). Paul Dupuis (Brown University) 26 October 2012
Examples One can use that π ( dx ) is the stationary distribution of the solution to √ dX = −∇ V ( X ) dt + 2 τ dW , as well as a variety of related discrete time models (will mention specific examples later). The function V ( x ) is defined on a large space, and includes, e.g., various inter-molecular potentials. In general, it may have a very complicated surface, with many deep and shallow local minima. A representative quantity of interest is � � V ( x ) e − V ( x ) /τ dx Z ( τ ) . Paul Dupuis (Brown University) 26 October 2012
Examples An example of a potential energy surface is the Lennard-Jones cluster of 38 atoms. This potential has ≈ 10 14 local minima. Paul Dupuis (Brown University) 26 October 2012
Examples An example of a potential energy surface is the Lennard-Jones cluster of 38 atoms. This potential has ≈ 10 14 local minima. The lowest 150 and their “connectivity” graph are as in the figure (taken from Doyle, Miller & Wales, JCP, 1999). Paul Dupuis (Brown University) 26 October 2012
Examples A second example. Other examples may have discrete (but very large!) domains of integration. A standard example is the Ising model (model for ferromagnetism). Paul Dupuis (Brown University) 26 October 2012
Examples A second example. Other examples may have discrete (but very large!) domains of integration. A standard example is the Ising model (model for ferromagnetism). Given a lattice (e.g., Λ = { 1 , . . . , N } × { 1 , . . . , N } ) and a neighborhood for each site i ∈ Λ (e.g., N ( i ) = { j : � j − i � 1 = 1 } ), define a measure on the collection of spins { + 1 , − 1 } Λ by e − β J ( x ) / Z ( β ) , where x i ∈ { + 1 , − 1 } and � � � J ( x ) = − J i , j x i x j − h i x i . i ∈ Λ j ∈ N ( i ) i ∈ Λ Paul Dupuis (Brown University) 26 October 2012
Examples A second example. Other examples may have discrete (but very large!) domains of integration. A standard example is the Ising model (model for ferromagnetism). Given a lattice (e.g., Λ = { 1 , . . . , N } × { 1 , . . . , N } ) and a neighborhood for each site i ∈ Λ (e.g., N ( i ) = { j : � j − i � 1 = 1 } ), define a measure on the collection of spins { + 1 , − 1 } Λ by e − β J ( x ) / Z ( β ) , where x i ∈ { + 1 , − 1 } and � � � J ( x ) = − J i , j x i x j − h i x i . i ∈ Λ j ∈ N ( i ) i ∈ Λ The analogue of the diffusion in the continuous state case are the Glauber dynamics (continuous time, finite state, jump Markov process). Paul Dupuis (Brown University) 26 October 2012
Monte Carlo approximation We focus on the continuous state problem. If X ( t ) is any discrete or continuous time Markov process with invariant measure µ , then under ergodicity � � � T T − 1 � µ T ( A ) = 1 or = 1 δ X ( t ) ( A ) dt δ X ( i ) ( A ) dt T T 0 i = 0 satisfies µ T ⇒ µ , and hence provides an numerical approximation. In practice one uses the discrete time process, and omits a “burn-in” period. Paul Dupuis (Brown University) 26 October 2012
Monte Carlo approximation Examples of discrete time processes. A simple example is the classical � Metropolis algorithm for π ( dx ) = e − V ( x ) /τ dx Z ( τ ) . It is based on two ideas: detailed balance and rejection . Paul Dupuis (Brown University) 26 October 2012
Monte Carlo approximation Examples of discrete time processes. A simple example is the classical � Metropolis algorithm for π ( dx ) = e − V ( x ) /τ dx Z ( τ ) . It is based on two ideas: detailed balance and rejection . For transition kernels with a density p ( x , dy ) = p ( y | x ) dy , detailed balance with respect to e − V ( x ) /τ dx means that e − V ( x ) /τ p ( y | x ) = e − V ( y ) /τ p ( x | y ) , and it automatically gives that e − V ( x ) /τ / Z ( τ ) is invariant with respect to p ( x , dy ) (integrate both sides w.r.t. x ). Paul Dupuis (Brown University) 26 October 2012
Monte Carlo approximation Standard Metropolis. Given the current location X ( i ) = x . Paul Dupuis (Brown University) 26 October 2012
Monte Carlo approximation Standard Metropolis. Given the current location X ( i ) = x . Choose a potential new location y according to a symmetric density 1 q ( y | x ) dy = g ( y − x ) dy (often uniform on the support of g ). Paul Dupuis (Brown University) 26 October 2012
Monte Carlo approximation Standard Metropolis. Given the current location X ( i ) = x . Choose a potential new location y according to a symmetric density 1 q ( y | x ) dy = g ( y − x ) dy (often uniform on the support of g ). Compute the acceptance probability 2 � 1 , e − V ( y ) /τ / e − V ( x ) /τ � � 1 , e − [ V ( y ) − V ( x )] /τ � A = min = min . Paul Dupuis (Brown University) 26 October 2012
Monte Carlo approximation Standard Metropolis. Given the current location X ( i ) = x . Choose a potential new location y according to a symmetric density 1 q ( y | x ) dy = g ( y − x ) dy (often uniform on the support of g ). Compute the acceptance probability 2 � 1 , e − V ( y ) /τ / e − V ( x ) /τ � � 1 , e − [ V ( y ) − V ( x )] /τ � A = min = min . Generate U ∼ Unif [ 0 , 1 ] , and accept the move if A ≥ U . 3 Paul Dupuis (Brown University) 26 October 2012
Monte Carlo approximation Standard Metropolis. Given the current location X ( i ) = x . Choose a potential new location y according to a symmetric density 1 q ( y | x ) dy = g ( y − x ) dy (often uniform on the support of g ). Compute the acceptance probability 2 � 1 , e − V ( y ) /τ / e − V ( x ) /τ � � 1 , e − [ V ( y ) − V ( x )] /τ � A = min = min . Generate U ∼ Unif [ 0 , 1 ] , and accept the move if A ≥ U . 3 Paul Dupuis (Brown University) 26 October 2012
Monte Carlo approximation Standard Metropolis. Given the current location X ( i ) = x . Choose a potential new location y according to a symmetric density 1 q ( y | x ) dy = g ( y − x ) dy (often uniform on the support of g ). Compute the acceptance probability 2 � 1 , e − V ( y ) /τ / e − V ( x ) /τ � � 1 , e − [ V ( y ) − V ( x )] /τ � A = min = min . Generate U ∼ Unif [ 0 , 1 ] , and accept the move if A ≥ U . 3 � 1 , e − [ V ( y ) − V ( x )] /τ � � e − V ( x ) /τ , e − V ( y ) /τ � Then e − V ( x ) /τ min = min plus symmetry of q guarantees that e − V ( x ) /τ / Z ( τ ) is invariant, but gives little information on rate of convergence. Paul Dupuis (Brown University) 26 October 2012
Recommend
More recommend
