Importance sampling for metastable and multiscale dynamical systems
Konstantinos Spiliopoulos
math.PR math.OC stat.ME






Importance sampling for metastable and multiscale dynamical systems
    K. Spiliopoulos[Department of  Mathematics and Statistics,
Boston University, Boston, MA, 02215, kspiliop@math.bu.edu. This work was partially supported by the National Science Foundation
CAREER award DMS 1550918 ]
    December 30, 2023
==========================================================================================================================================================================================================================




In this article, we address the issues that come up in the design of importance sampling schemes for rare events associated to stochastic dynamical systems. We focus on the issue of metastability and on the effect of multiple scales.  We discuss why seemingly reasonable schemes that follow large deviations optimal paths may perform poorly in practice, even though they are asymptotically optimal. Pre-asymptotic optimality is important when one deals with metastable dynamics  and we discuss possible ways as to how to address this issue. Moreover, we discuss how the effect of the multiple scales (either in periodic or random environments) on the efficient design of importance sampling should be addressed. We discuss the mathematical and practical issues that come up, how to overcome some of the issues and discuss future challenges.





§ INTRODUCTION




In this paper, we discuss recent developments on importance sampling methods for metastable dynamics that may also have multiple scales. Development of accelerated Monte Carlo methods for metastable, multiple-scale processes is of great interest.  Importance sampling is a variance reduction technique in Monte-Carlo simulation, which is especially relevant when dealing with rare events. Since its introduction, importance sampling has been one of the most popular techniques for rare event simulation. There is a vast literature of papers investigating its applications from a broad set of sciences including engineering, chemistry, physics, biology, finance, insurance, e.g., <cit.>.

Consider a sequence {X^ϵ}_ϵ>0 of random elements and assume that we want to estimate the probability 0<p^ϵ=ℙ[X^ϵ∉𝒟∪∂𝒟]≪ 1 for a given set 𝒟, such that the event {X^ϵ∉𝒟∪∂𝒟} is unlikely for small ϵ. If closed form formulas are not available, or numerical approximations are either too crude or unavailable, then one has to resort in simulation. It is well known that standard Monte-Carlo simulation techniques (i.e., using the unbiased estimator p̂^ϵ=1/N∑_j=1^N1_X^ϵ,j∉𝒟∪∂𝒟) perform rather poorly in the rare-event regime. As it is known, see for example <cit.>,  in order to  achieve relative error smaller than one using standard Monte Carlo, one needs an effective sample size N≈ 1/p^ϵ. In other words, for a fixed computational
cost, relative errors grow rapidly as the event becomes more rare.   Thus  standard Monte-Carlo is infeasible for rare-event simulation.




The goal of importance sampling is to simulate the system under an
alternative probability distribution ℙ̅ instead of the original probability ℙ. Let's say for example that we are interested in the estimation of
 
    𝔼_y[e^-1/ϵh(X^ϵ_T)]
    or ℙ_y[τ^ϵ_𝒟∪∂𝒟≤ T]

where h:ℝ^d↦ℝ is a positive function, T>0, ϵ>0, y∈ D is the initial point, τ^ϵ_𝒟∪∂𝒟 is exit time from the set 𝒟∪∂𝒟, X^ϵ is a stochastic process modeling the dynamics.  Also, notice that the probability  above can be considered (modulo the important technical point of lack of continuity) as a special case of 𝔼_y[e^-1/ϵh(X^ϵ_T)], when h is for example chosen such that h(x)=0 for x∉𝒟∪∂𝒟 and h(x)=+∞ for x∈𝒟∪∂𝒟.

When rare events dominate, then standard Monte-Carlo
methods perform poorly in the small noise limit. 






 Then, to estimate  𝔼_y[e^-1/ϵh(X^ϵ_T)], one generates iid samples X^ϵ_(k) from ℙ̅ and
uses the importance sampling estimator 

    1/N∑_k=1^N  e^-1/ϵh(X^ϵ_(k))dℙ/dℙ̅(X^ϵ_(k)).

The key question is the design of ℙ̅ such that the second moment 𝔼̅_y[e^-1/ϵh(X^ϵ_T)(dℙ/dℙ̅)(X^ϵ_·)]^2




(and hence the variance) is minimized. 𝔼̅ is
the expectation operator under ℙ̅. The choice of the appropriate alternative measure ℙ̅  is closely related to certain Hamilton-Jacobi-Bellman (HJB) equations.


The first issue that we address is the effect of rest points (and metastability in general) on importance sampling. It turns out that when dealing with metastability, even seemingly reasonable schemes that are also asymptotically optimal, may perform poorly in practice. This includes also changes of measure that try to enforce the simulated trajectories to follow large deviations most likely paths. The reason for the degradation in performance is the role of prefactors. Prefactors can become very important when rest points are included in the domain of interest for the simulation. Large deviations based change of measures may not account for the prefactors, as they rely on logarithmic asymptotics. We elaborate on these issues and discuss potential ways on how the issue can be addressed.

The second issue that we address is the effect of multiple scales on the design of provably-efficient importance sampling methods. It turns out that when the dynamical system has widely separated multiple scales, then one can use averaging and homogenization techniques. However, as we shall see, it is not sufficient to base the design of importance sampling on the effective homogenized dynamics. The local information needs to be taken into account. Mathematically this is done using the so called cell problem, or macroscopic problem, in the theory of periodic and random homogenization.

The rest of the article is summarized as follows. In Section <ref> we review the classical large deviations theory and the setup of importance sampling for small noise diffusions. In Section <ref> we discuss the effects of rest points, i.e. of stable and unstable equilibrium points, in the design of importance sampling. We argue why asymptotic optimality may actually not mean good practical performance and we also argue that following large deviations most likely optimal paths may lead to poor performance. In addition, we present constructions that lead to guaranteed good performance. We supplement the theoretical arguments by simulation studies. We refer the interested reader to <cit.> for more details. In Sections <ref> and <ref>, we address the design of importance sampling schemes in the presence of multiple scales. We construct asymptotically optimal schemes in the presence of multiple scales. To be more precise, in Section <ref> we consider overdamped Langevin dynamics in periodic multiscale environments and we review the related large deviations theory and importance sampling theory, presenting simulation studies.  The interested reader can also consult <cit.>. In Section <ref> we review recent developments in large deviations and importance sampling for multiscale dynamics in random environments, see also <cit.>. In Section <ref> we describe how one can combine the results of Section
<ref> with those of Sections <ref> and <ref> and also review future directions.

For the sake of concreteness and for exposition purposes we restrict the presentation of this article in the case of diffusions with gradient drift and constant diffusivity, which also implies reversible diffusion dynamics. However, we mention that almost all of the arguments can and have been generalized to the case with general state dependent drift and diffusion coefficient, especially those about the effect of multiple scales on importance sampling, see <cit.>. For results in the infinitely dimensional case we refer the interested reader to <cit.>.




§ REVIEW OF LARGE DEVIATIONS AND IMPORTANCE SAMPLING THEORY FOR DIFFUSIONS

Let us briefly review the setup for small noise diffusions in ℝ^d (e.g. <cit.>) without the effect of multiple scales. Let W_t be a standard d-dimensional Wiener process and consider

    dX_t^ϵ=-∇ V( X_t^ϵ)  dt+√(ϵ)Γ dW_t,X_t_0^ϵ=y.



Large deviations principle for the process X^ϵ_t is well known (e.g, <cit.>). In particular, the action functional  for the process
X^ϵ_t, t_0≤ t≤ T, in 𝒞([t_0,T]) as ϵ↓ 0 has the form 1/ϵS_t_0T(ϕ), where 






    S_t_0T(ϕ)=1/2∫_t_0^T(ϕ̇_s+∇ V(ϕ_s))^T[ΓΓ^T]^-1(ϕ̇_s+∇ V(ϕ_s))ds,    if ϕ∈𝒜𝒞([t_0,T])
     
    
                   +∞,    otherwise.

Here 𝒞([t_0,T]), 𝒜𝒞([t_0,T]) are the sets of continuous and absolutely continuous functions on [t_0,T] respectively.

Then, under fairly general conditions,

    𝔼_y[e^-1/ϵh(X^ϵ_T)]≈ e^-1/ϵinf{S_t_0T(ϕ)+h(ϕ_T):ϕ,ϕ_t_0=y},  as ϵ↓ 0.


A simple application of Jensen's inequality together with Varadhan's integral lemma (e.g., <cit.>) shows that an asymptotically optimal ℙ̅ should satisfy

    lim_ϵ→ 0ϵln𝔼̅[e^-1/ϵh(X^ϵ_T)dℙ/ dℙ̅]^2=-2G(t_0,y),  with  G(t,x)=inf_ϕ∈𝒜𝒞([t,T]), ϕ_t=x{S_tT(ϕ)+h(ϕ_T)}

Turning to importance sampling, for ℙ̅ that are absolutely continuous with respect to ℙ, Girsanov's formula implies 


    dℙ̅/dℙ=e^-1/2ϵ∫_0^T|v_s|^2
    ds+1/√(ϵ)∫_0^Tv_sdW_s

where v_t is a progressively measurable
process (control) such that the right hand side  is a martingale (with respect to an appropriate filtration). Under ℙ̅, X^ϵ satisfies 

    dX_t^ϵ=[-∇ V( X_t^ϵ)+Γ v_t]  dt+√(ϵ)Γ dW̅_t,  with  W̅_t=W_t-1/√(ϵ)∫_t_0^tv_ρdρ


So, the problem is restricted to choosing the control v_t  optimally (i.e., such that the second moment is minimized)  and then using the estimator
based on iid samples generated from ℙ̅ under (<ref>).






Under appropriate conditions, the zero-variance (i.e. the best) change of measure is based on the control v_t given by the formula v_t=u̅(t,X^ϵ_t) where

v̅(t,x)=-Γ^T∇ G^ϵ(t,x) 

where G^ϵ(t,x), with terminal condition G^ϵ(T,x)=h(x),  is the  solution to the PDE, of HJB type:

    ∂_tG^ϵ(t,x)-∇ V(x)·∇ G^ϵ(t,x)-1/2|Γ^T∇ G^ϵ(t,x)|^2+ϵ/2tr[ΓΓ^T∇^2G^ϵ(t,x)]=0.


Since (<ref>) is not tractable, it is standard approach to go to the viscosity limit ϵ↓ 0. Then G(t,x)=lim_ϵ↓ 0G^ϵ(t,x)  is the viscosity solution to the  HJB equation with Hamiltonian

    H(x,p)=<-∇ V(x), p>-1/2‖Γ^T p‖^2

 i.e., to the equation
     
    ∂_tG(t,x)-∇ V(x)· DG(t,x)-1/2|Γ^T DG(t,x)|^2=0,  G(T,x)=h(x).


Notice that by control arguments, e.g., see <cit.>, we can also write

    G(t,x)=lim_ϵ↓ 0G^ϵ(t,x)=inf_ϕ∈𝒜𝒞([t,T]), ϕ_t=x{S_tT(ϕ)+h(ϕ_T)}.


In fact, more is true. A smooth function U̅(t,x):[0,T]×ℝ^d↦ℝ is called a subsolution to the HJB
equation (<ref>) with ϵ=0 if

    ∂_tU̅(t,x)-∇ V(x)·∇U̅(t,x)-1/2|Γ^T∇U̅(t,x)|^2≥ 0, U̅(T,x)≤ h(x).


It turns out (Theorem 4.1 in <cit.>), that appropriate, smooth subsolutions are enough.
If U̅(t,x)∈𝒞^1,1([t_0,T]×ℝ^d) satisfies (<ref>) and
the feedback control to use in (<ref>) is v_t=-Γ^T∇U̅(t,X^ϵ_t), then




    G(t_0,y)+U̅(t_0,y)≤lim inf_ϵ→0-ϵln𝔼̅[e^-1/ϵh(X^ϵ_T)dℙ/dℙ̅]^2≤ 2G(t_0,y).



Therefore,  asymptotic optimality is attained if U̅ satisfies U̅(t_0,y)=G(t_0,y)=lim_ϵ↓ 0G^ϵ(t_0,y) since then lower and upper bound agree. The design and analysis of importance sampling schemes based on the systematic connection with subsolutions to the appropriate HJB and Isaacs  equations goes back to <cit.>. See also  <cit.> for the closely related concept of Lyapunov inequalities.



The importance sampling simulation scheme in order to estimate θ^ϵ(t_0,y)≐𝔼_t_0,y[  e^-1/ϵh(X_T^ϵ
)] goes as follows. Let X^ϵ,v be
the solution to the SDE

    dX_t^ϵ,v   =(-∇ V( X_t^ϵ,v)+Γ v_t) dt+√(ϵ)Γ dW_t,X_t_0^ϵ,u=y.




  * Consider v_t=u̅(t,X_t^ϵ,v)=-Γ^T∇_xU̅(t,X_t^ϵ,v)  with U̅ an appropriate subsolution, i.e., it satisfies (<ref>)

  * Consider the estimator

    θ̂^ϵ(y)≐1/N∑_j=1^N[  e^-1/ϵh(X_T^ϵ,v(j))Z_j^v]

where

    Z_j^v≐ e^-1/2ϵ∫_0^Tu̅(
    t,X_t^ϵ,v(j))   ^2dt-1/√(ϵ)
    ∫_0^T<u̅(  t,X_t^ϵ,v(j)),  dW_t
    (j)>

and (W(j),X^ϵ,v(j)) is an independent sample generated from
(<ref>) with control v_t=u̅(  t,X_t^ϵ
,v(j)).



We conclude this section, with the remark that a choice of the control v_t based on a subsolution as defined by (<ref>)  only guarantees logarithmic asymptotic optimality and does not say something about the important effect of pre-factors. As we will see in Section <ref>, this can imply degradation in the performance of the algorithm in problems with metastability.  When dealing with metastability issues, things may be even more problematic if one is using the exact solution to the association HJB equation, G(t,x). While this may be not be a problem for problems that do not involve rest points (i.e. does not involve stable or unstable equilibrium points) in the domain of interest, it does become problematic when dealing with metastability issues.


Obtaining  accurately the solution G(t,x) to the HJB equation (<ref>), analytical or numerical,  is challenging in high dimensions. However, even if this were possible, the solution by itself is not always  suitable for importance sampling when one is interested in computing escape or transition probabilities.
The issue is that in these cases, the solution is a viscosity solution with a discontinuous derivative at the rest point (stable or unstable equilibrium points) and with negative definite generalized second derivative there. Physically, the exact solution to the HJB equation attempts at each point in time and space to force the simulated trajectories to follow a most likely large deviations optimal path. However, by standard control arguments, see <cit.>, the discontinuity of the spatial derivative at the rest point, implies that multiple most likely optimal paths exist. As a consequence, the noise
can cause trajectories to return to a neighborhood of the origin, thereby
producing large likelihood ratios.  In Section <ref>, we will see that this is a serious issue, leading to poor performance, even in dimension one where one can solve the HJB equation analytically. Importance sampling, when dealing with state dependent metastable dynamical systems, needs to be addressed from a global point of view and not local.







§ THE EFFECT OF REST POINTS  ON IMPORTANCE SAMPLING







As it is shown, mathematically and numerically, in <cit.>, in dynamical systems that exhibit metastable behavior standard simulation methods do not readily apply. Asymptotic optimality is necessary but not sufficient for good performance
due to the non-trivial effect of the pre-factors.  The pre-factor computations in <cit.> prove that there is non-trivial interaction of parameters such as
the strength of the noise ϵ and the terminal time T. We remark here that this is in contrast to
escape probabilities for other well studied problems, such as stochastic networks, e.g., <cit.>,
because there the proximity of the rest point has little impact on either the asymptotic rate of decay or the pre-exponential term.

These interactive effects vanish in the logarithmic limit as the noise goes to zero, but they have
a significant effect on the performance of the algorithms. The following question immediately presents itself:


  * Is it sufficient to have schemes that are only asymptotically logarithmical optimal, in the sense that the second moment of the estimator satisfies (<ref>)? What about pre-factors? Are they truly negligible in practice in the rare event regime?

  * Can we construct a subsolution U̅(t,x) that not only satisfies (<ref>) but it also takes care of the prefactor effects?




 §.§ Effects in the prelimit

Let us demonstrate the effect of prefactors on the behavior of estimators in the following classical simple setting. Let us assume that the diffusion coefficient Γ=I, and that x=O is the global minimum for V(x). In particular, let us assume that DV(O)=0 and that DV(x)≠ 0 for every x≠ O. Define

    𝒟={ x∈ℝ^d: 0≤ V(x)<L }

and let A_c={ x∈ℝ^d: V(x)=c}. Then for an initial point y such that 0≤ V(y)<L, let us assume that we want to estimate

    θ^ϵ(t,y)=ℙ_t,y{ X^ϵ hits  A_L before time T}.


A classical quantity if interest in metastability theory is the quasipotential, see <cit.>. The quasipotential with respect to the equilibrium point O is defined as follows

    W(O,x)={ S_0T(ϕ): ϕ∈𝒞([0,T]), ϕ(O)=0, ϕ(T)=x, T∈(0,∞)}


Under the assumptions of this section, the  quasi-potential is computable in closed form <cit.>:

W (O,x) = 2V(x)  for  x∈{y∈𝒟∩∂𝒟: V(y)≤inf_z∈∂𝒟V(z)}.






Now, if we define τ^ϵ=inf{t>0: X_t^ϵ∉𝒟}, then, as it is shown in <cit.> we have that

 lim_ϵ↓ 0ϵln𝔼τ^ϵ=inf_z∈∂𝒟W(O,z).

Thus, the quasi-potential allows to approximate exit times in the logarithmic large deviations regime,  <cit.>.  Many quantities in the theory of metastability are defined via the quasi-potential.
The quasi-potential characterizes the leading asymptotics of exit times and exit probabilities, approximates transition rates for reversible and irreversible systems and allows to qualitatively describe transitions between stable attractors if the system has many of them; see also <cit.> for more details. These conclusions hold for both gradient and non-gradient cases, but in the gradient case the quasi-potential is computable in closed form.

Turning now to importance sampling,  it is easy to verify that the quasi-potential is a stationary subsolution to the associated HJB equation (<ref>) with ϵ=0, by adding an appropriate constant C in order to justify the necessary boundary and terminal conditions. In particular, U̅_QP(x)=2L-W(O,x) defines a subsolution for (<ref>). It turns out, see <cit.>, that the quasipotential yields a  reasonable change of measure if rest points are not part of the domain of interest. However, this is no longer true if rest points are included in the domain of interest.

Let us denote Q^ϵ(0,y;u̅)=𝔼̅[e^-1/ϵh(X^ϵ_T)dℙ/ dℙ̅]^2 to be the second moment of the estimator constructed using the control u̅. Based now on the arguments of <cit.> one can prove the following representation for the second moment of the estimator estimator based on the change of measure induced by the control u̅(t,x)=-∇U̅_QP(x)

    -ϵlog Q^ϵ(0,y;u̅)=inf_v∈𝒜
    𝔼[  1/2∫_0^τ̂^ϵ‖
    v(s)‖ ^2ds-∫_0^τ̂^ϵ‖u̅(X̂^ϵ_s)‖ ^2ds+∞1_{τ̂^ϵ>T}]  .

 where X̂^ϵ_s is the unique solution to the SDE
 
    dX̂^ϵ_s=-DV(X̂^ϵ_s)ds+[
    √(ϵ)dW_s-[u̅(X̂^ϵ_s)-v(s)]ds]

with initial condition X̂^ϵ_0=y and τ̂^ϵ is the first time that X̂^ϵ exits from 𝒟.

It is important to note that (<ref>) provides a non-asymptotic representation for the second moment of the estimator. By the arguments of <cit.>, we can choose a particular admissible control v(s) in (<ref>) so that the following takes place. Let T be large and let 0<K<T so that the time interval [0,T] is split into [0,T-K) and [T-K,T]. Set v(s)=0 for s∈[0,T-K). The resulting dynamics for X̂^ϵ is stable for s∈[T-K,T] and with high probability the process will stay around the point y for s∈[0,T-K). In the time interval [T-K,T], we set v(s) so that escape happens prior to T. Then, it can be shown that there are positive constants C_1,C_2<∞, so that

    Q^ϵ(0,y;u̅)≥ e^-1/ϵC_1+C_2(T-K).


This bound indicates that if T is large, one may need to go to considerably small values of ϵ in order to achieve the theoretical optimal asymptotic performance. We also remark that if T is large (see Chapter 4 of <cit.>), G(0,y) and U̅(y) get closer in value. Thus, by (<ref>) and for large enough T, the particular importance sampling scheme is asymptotically optimal.

Hence, we have just seen an example where an importance sampling estimator is almost asymptotically optimal, but it does not perform that well pre-asymptotically due to the effect of the possibly long time horizon T and its interplay with ϵ.



 §.§ The problems arising when following large deviations asymptotically most likely paths and a remedy to the problem

The connection of change of measures with HJB equations via large deviations is well situated for a systematic treatment of dynamic importance sampling schemes  for state dependent processes like diffusions (<ref>). For small noise diffusions the theoretical framework of subsolutions to HJB equations and their use for Monte Carlo methods can be found in <cit.>. It was a common belief for sometime that if the underlying stochastic process has a large deviations principle and if the change of measure is consistent with
 the large deviations asymptotically most likely path leading to the rare event (an open-loop control), then the resulting importance sampling scheme would be optimal. 

However, such heuristics have been shown to be unreliable in general and simple examples have been
constructed showing the failure of the corresponding importance schemes even in very simple settings <cit.>. This is due to the presence
of “rogue-trajectories", i.e., unlikely trajectories, that are likely enough to increase likelihood ratios to the point that
 the performance is comparable to standard Monte Carlo. This is especially true for metastability problems (i.e., when transitions between fixed points occur at suitable (large) timescales)
 where multiple nearly optimal paths may exist.

 Use of dynamic changes of measure, i.e. based on feedback controls (time and location dependent) becomes important, see <cit.>. However, even changes of measures that are based on feedback controls, that are consistent with large deviations and lead to asymptotically optimal change of measures can also be problematic in practice. We demonstrate this below in Table <ref>. Namely, as it turns out, in the presence of rest points and metastability, the prefactors may affect negatively the behavior of estimators even if one is using asymptotically optimal changes of measure in the spirit of (<ref>). Hence, it becomes important to use dynamic change of measures that are based on subsolutions but lead to good performance even pre-asymptotically.


To that end, novel explicit simulation schemes are then constructed in <cit.> that perform provably-well
both asymptotically and non-asymptotically, even when the  simulation time is long.
These constructions are based on large deviations asymptotics <cit.>, stochastic control arguments and asymptotic expansions <cit.>  and detailed asymptotic analysis of the subsolution to the
associated HJB in the neighborhood of the rest point where the potential can be thought of as being approximately quadratic. Essentially, due to the fact that near the rest point, the potential can be thought of as being approximately quadratic, one can hope to solve or to approximate the solution to the associated variational problem there. Then one needs to patch this solution together with the quasipotential based subsolution (which is a good subsolution away from the rest point) in the right way. Then, the combined subsolution, see U̅^δ(t,x) in (<ref>), turns out to be a good approximation to the zero variance change of measure. Such schemes lead to importance sampling algorithms with provably-good performance for all small ϵ>0 and without suffering from bad prefactor effects.

In order to illustrate the point, let us briefly demonstrate such a construction in the case of dimension one, see <cit.>. So, let us assume that V(x)=λ/2x^2 with λ>0 and let us assume that we study the problem of crossing a level set, say L, of the potential function V(x). Here, we can compute G(t,x) in closed form and we get

    G(t,x)      =inf_ϕ_t=x,V(ϕ_T)=L{1/2∫_t
    ^T‖ϕ̇_s+λϕ_s‖ ^2ds}
      =inf_x̂∈ V^-1(L)λ(  x̂-xe^λ
    (t-T))  ^2/1-e^2λ(t-T).


Notice, that G(t,x) is also a viscosity solution to the ϵ=0 HJB equation  (<ref>) when supplemented with the appropriate boundary conditions. Hence, based on (<ref>) a change of measure based on G(t,x), i.e., using the control u(t,x)=-∂_xG(t,x), is expected to yield an asymptotically efficient estimator. While this is true, we will see below that this is not sufficient to yield good performance. The fact that the function G(t,x) is not continuously differentiable in the domain of interest, implies that multiple optimal paths exist, which is an intuitive reason for the degradation in performance that will be demonstrated below.

However, by appropriately mollifying  G(t,x) and combining it with the quasipotential subsolution (as constructed in Section <ref>), one can construct a global subsolution which performs provably well even pre-asymptotically. The point  is that G(t,x) provides a good change of measure
while near the rest point, whereas the quasipotential induced subsolution U̅_QP(x)=2L-W(O,x) provides a
good change of measure away from the rest point. There are a few more issues
to deal with though. The first one is that G(t,x) is discontinuous near
t=T. The second one is that we need to put them together in a smooth way
that will define a global subsolution.

Since G(t,x) is discontinuous at t=T,  we introduce two
mollification parameters t^∗ and M that will be appropriately chosen
as functions of ϵ. Motivated by the fact that G(t,x) is a good subsolution near the equilibrium point, we fix another parameter
L̂∈(0,L]. In the one-dimensional case, it is easy to solve the
equation V(x^*)=L̂ and in particular we get that x^*=±x̂ where x̂
=√(2L̂/λ). As a matter of fact, instead of using G(t,x) directly,
we set

    F^M(t,x;x̂)       = λ(  x̂-xe^λ
    (t-T))  ^2/1/M+1-e^2λ(t-T)



In order now to pass smoothly between the U̅_QP(x) and F^M(t,x;x̂) or F^M(t,x;-x̂) without violating the subsolution property, we use the
exponential mollification, see <cit.>

    U^δ(t,x)=-δlog(  e^-1/δU̅_QP(x)+e^-1/δ[  F^M(t,x;x̂)+U̅_QP(x̂)]  +e^-1/δ[  F^M(t,x;-x̂)+U̅_QP(-x̂)]  )

It is easy to see that as δ↓ 0

    lim_δ↓ 0U^δ(t,x)=min{U̅_QP(x), F^M(t,x;x̂), F^M(t,x;-x̂)}


Clearly, if we choose L̂=L, then we get U̅_QP(x̂)=0. Based on these constructions, a provably efficient importance sampling scheme is constructed in
<cit.>, based on the
subsolution

    U̅^δ(t,x)={[ U̅_QP(x),   t>T-t^∗; U^δ(t,x),  t≤ T-t^∗ ].


It turns out that U̅^δ(t,x) is a global smooth subsolution which
has provably good performance both pre-asymptotically and asymptotically.  The role of the exponential mollification is to allow a smooth transition between the region that is near the equilibrium point and the region that is far away from it. 




















































































The
precise optimality bound and its proof guide the choice of the parameters
δ,t^∗,M and L̂.
For the convenience of the reader, we present in Table <ref> the suggested values for (δ,L̂,M,t^*), given the value of
the strength of the noise ϵ>0.

 

We refer the interested reader to <cit.> for further details on the theoretical performance of the algorithm and on the choice of parameters.

In order to illustrate in a simple setting the effect of prefactors in the presence of metastable effects, we record in Table <ref> Monte Carlo estimates based on K=10^7 trajectories for the exit time distribution ℙ_y[τ^ϵ_𝒟∪∂𝒟≤ T]  from the basin of attraction of the left attractor of
the potential of Figure 1 for the process X^ϵ given by (<ref>) with Γ=I. We used the importance sampling (IS) methods of <cit.>, i.e., the change of measure based on the subsolution (<ref>) and record estimates for different pairs (ϵ,T).


In the figures next to Table <ref>, we  compare the relative errors per sample of (a): the algorithm, which is optimal for all ϵ>0, i.e the one based on the subsolution U̅^δ(t,x), with (b): the IS algorithm that is consistent with the
large deviations asymptotically most likely path leading to the rare event, i.e the one based on the actual solution G(t,x) of the associated HJB equation. Notice  however that the IS algorithm based on G(t,x) is only asymptotically optimal in the large deviations logarithmic sense as ϵ↓ 0 (i.e., it satisfies (<ref>)).








































Using relative error per sample as comparison criterium, we compare the two algorithms for two values of ϵ, one for which the events are not so rare (ϵ=0.13) and one for which the events are very rare (ϵ=0.05). Exact values are in the table, and we remark for completeness that intermediate behavior is qualitatively the same. Both algorithms perform well when T is small, but the algorithm that is based on the solution of the associated HJB equation, which is only logarithmic asymptotically optimal,
starts deteriorating considerably as T gets large. The latter  is an effect of the pre-factors becoming important. On the other hand, the change of measure
constructed in <cit.> that takes into account the pre-factor information and is pre-asymptotically optimal, yields optimal performance independently of the values ϵ and T with relative errors around one, meaning that the values recorded at the table are reliable. It is important to note that due to large deviations, exit happens in long time scales, which implies that reliable estimates, especially when T is large, are essential.




§ IMPORTANCE SAMPLING FOR ROUGH ENERGY LANDSCAPES

In Section <ref>, we reviewed some of the practical issues that come up when  one is trying to apply importance sampling techniques to metastable dynamics. While in Section <ref> we ignored the effect of multiple scales, the goal of this section is to address the role of multiple scales in the design of asymptotically optimal importance sampling schemes.

 A particular model of interest in chemical physics is the first order Langevin equation (<ref>).  Let us consider

    dX_t^ϵ,δ=[ -ϵ/δ∇ Q(
    X_t^ϵ,δ/δ) -∇ V( X_t^ϵ,δ) 
    ] dt+√(ϵ)√(2D)dW_t,X_0^ϵ,δ=y, , 0<ϵ,δ≪ 1,

where  the two-scale potential is composed
by a large-scale part, V(x), and a fluctuating part, ϵ Q(x/δ). If Q is periodic then we have a periodic environment, whereas if Q is random then we have a random environment. Models like (<ref>) can be used to model rough energy landscapes
<cit.>.  As it has been suggested
(e.g., <cit.>), the associated energy landscapes of certain
biomolecules can be rugged (i.e., consist of many local “small" minima within local deep minima separated
by barriers of varying heights). When one is interested in rare events, large deviations and Monte Carlo methods are relevant.




If Q(y) is periodic,  large deviations for multiscale diffusions in periodic environments are obtained in <cit.> for all possible interactions between ϵ and δ, setting the ground for the mathematical formulation of the related importance sampling theory, <cit.>. The novel feature is that the optimal change of measure for importance sampling is not based only on the gradient of the homogenized HJB equation (as in Subsection <ref>). The effect of fluctuations, which is quantified via the solution to the “cell problem" in  homogenization <cit.>,
is equally important. The cell problem is the solution to a Poisson type PDE. It is used to define the so called “corrector", which characterizes the first order correction in the approximation of the multiscale HJB  by its homogenized limit.  Therefore, when compared to the case without multiple scales, one needs more detailed information in order to guarantee, at least, asymptotic optimality.



For example, consider model (<ref>) in the case ϵ/δ↑∞. Define the Gibbs measure

    μ(dy)=1/L e^-Q(y)/Ddy,   L=∫_𝕋^d e^-Q(y)/Ddy.


Then denote by χ(y)  the smooth solution to the “cell problem”

    -∇ Q(y)·∇χ(y)+ Dtr[∇^2χ(y)] = ∇ Q(y),    ∫χ(y)μ(dy) = 0.



The following large deviations result holds which is a special case of the results of <cit.>. In particular, <cit.> covers the case of general state dependent drift (not necessarily of gradient form) and state dependent diffusion coefficient.

 Assume that the functions ∇ Q(
y ) and ∇ V(x) are continuous and globally bounded,
as are their partial derivatives up to order 1 in y and order 2 in x respectively. Let {X^ϵ,δ,ϵ,δ>0} be the unique
strong solution to (<ref>). Let

    r(x)       =-∫_𝕋^d(  I+∂χ(y)/∂
    y)  μ(dy) ∇ V(x),
    
    q       =2D ∫_𝕋^d(  I+∂χ(y)/∂
    y)  (  I+∂χ(y)
    /∂ y)^Tμ(dy),

where I denotes the identity matrix. If ϵ/δ→∞, then {X^ϵ,δ,ϵ,δ>0} converges in probability as ϵ,δ→ 0 to the solution of the ODE

    dX̅_t=r(X̅_t)dt

and satisfies a large deviations principle with rate function

    S_tT(ϕ)=
    1/2∫_t^T(  ϕ̇_s-r(ϕ_s)) q^-1(  ϕ̇_s-r(ϕ_s))^T ds    if ϕ∈𝒜𝒞
    ([t,T]),ϕ_t=x
    
    +∞   otherwise.




In addition, it turns out that an asymptotically efficient change of simulation measure can be constructed analogously to Section <ref>, but  based on the  feedback control (see Theorem 4.1 in <cit.>)

    v_t=u̅(t,X^ϵ_t,X^ϵ_t/δ),  with  u̅(t,x,y)=-√(2D)( I+ ∂χ(y)/∂ y) ^T∇_xU̅(t,x).


U̅(t,x) satisfies the inequalities in (<ref>) with the homogenized (averaged) coefficients r(x) and q in place of the original ones -∇ V(x) and Γ= √(2D) I (compare with (<ref>)). In particular, the second moment of an estimator with change of measure based on the control v_t by (<ref>) will satisfy (<ref>); this is Theorem 4.1 in <cit.>.



Thus, compared to the case without
multiscale features, one needs to consider the extra factor ( I+ ∂χ(y)/∂ y), that can be thought as the appropriate weight function, to achieve asymptotic optimality. In the absence of multiple scales, i.e., when Q=0, we have χ=0 and we recover the
case studied in Section <ref>. The numerical simulation studies of <cit.> verify the need for accounting for the local environment via the weights ( I+ ∂χ(y)/∂ y) in the change of simulation measure.

Before illustrating the performance of this importance sampling scheme in a simulation study, let us demonstrate theoretically the necessity to include the cell problem information in the design of the change of measure. For simplicity purposes, let us restrict attention to dimension one. As we have seen before, the effective diffusion coefficient is given by

    q=2D ∫_𝕋(  1+∂χ/∂ y) ^2μ(dy)


In this case, the optimal change of measure is based on the control

    u̅(t,x,y)=-√(2D)( 1+ ∂χ(y)/∂ y) ∂_xU̅(t,x).

 So, let us assume that one is using instead the change of measure, based on the control dictated by the averaged dynamics. Namely, let us assume that the control in question is
û(t,x)=-√(q)∂_xU̅(t,x).

A verification theorem, see <cit.> for details,  shows that one would need a statement of the form

    "𝔼∫_t^T[  √(2D)(  1+∂χ/∂ y( X^ϵ,δ_s/δ))  -√(q)]
    ds→ 0"


By averaging principle, this is true if

    √(q)=∫_𝕋√(2D)(  1+∂χ(y)/∂ y) μ(dy).


However, this is impossible, since

    (∫(  1+∂χ(y)/∂ y)  μ(dy))^2≠∫(  1+∂χ(y)/∂ y)^2μ(dy).


This last property explains mathematically why, the local information, as quantified via the cell problem, needs to be taken into account in the design of importance sampling. In Section <ref>, we will also see numerical evidence of this issue.


 §.§ A simulation study


Let us demonstrate the performance of the importance sampling scheme in a simple simulation study. Consider the one well potential function with diffusion coefficient D=1,

    V(x)=1/2x^2, 
    Q(y)=cos(y)+sin(y)





Assume that we want to estimate θ(ϵ,δ)=𝔼[e^-1/ϵh(X^ϵ,δ_1)], where

 h(x)=(|x|-1)^2.


It is easy to see that we are dealing with a rare event here, as the function h(x) is minimized at |x|=1. Let us compare the following three different estimators

    θ̂_0(ϵ,δ)   =   1/K∑_j=1^K[e^-1/ϵh(X^ϵ,δ_1(j))]  - standard Monte Carlo 
    θ̂_1(ϵ,δ)   =   1/K∑_j=1^K[e^-1/ϵh(X̅^ϵ,δ, u̅_1(j))Z^u̅_j]  - optimal
    θ̂_2(ϵ,δ)   =   1/K∑_j=1^K[e^-1/ϵh(X̅^ϵ,δ, û_1(j))Z^û_j]  - ignores local information


where we have defined the controls


  * u̅(t,x,y)=-√(2)( 1+ ∂χ(y)/∂ y)G_x(t,x)–asymptotically optimal.

  * û(t,x)=-√(q)G_x(t,x)–based only on the homogenized system.

and the likelihood ratio is Z^u_j=dP/dP̅(X̅^ϵ,δ, u_1(j)). Notice that in this case, we can compute

    1+∂χ(y)/∂ y=e^Q(y)/∫_𝕋 e^Q(y)dy,

which justifies the interpretation of the term 1+∂χ(y)/∂ y as the proper weight term needed that takes into account the local information.

In Table <ref>, we see simulation studies based on N=10^7 simulation trajectories each, for the estimation of θ(ϵ,δ) using the three different estimators. The measure of comparison is chosen to be the relative error per sample, defined to be

    ρ̂_i(ϵ,δ)≐√(N)√(Var(θ̂_i(ϵ,δ)))/θ̂_1(ϵ,δ).





It is clear, that the importance sampling scheme based on the asymptotically optimal change of measure u̅(t,x,y) outperforms the standard Monte Carlo estimator in which no change of measure is being done. It also outperforms,  the estimator based solely on the homogenized system, which ignores the local information characterized by solution to the cell  problem χ(y).




§ IMPORTANCE SAMPLING FOR MULTISCALE DIFFUSIONS IN RANDOM ENVIRONMENTS


Let 0<ϵ,δ≪ 1 and consider the process (X^ϵ, Y^ϵ)={(X^ϵ_t, Y^ϵ_t), t∈[0,T]} taking values in the space ℝ^m×ℝ^d-m that satisfies the system of SDEs


    dX^ϵ_t   =   [  ϵ/δb(Y^ϵ_t,γ)+c(  X^ϵ_t
    ,Y^ϵ_t,γ)]   dt+√(ϵ)
    σ(  X^ϵ_t,Y^ϵ_t,γ)
    dW_t,
    
    dY^ϵ_t   =   1/δ[  ϵ/δf(Y^ϵ_t,γ)  +g(  X^ϵ_t
    ,Y^ϵ_t,γ)] dt+√(ϵ)/δ[
    τ_1(  Y^ϵ_t,γ)
    dW_t+τ_2(Y^ϵ_t,γ)dB_t], 
    
    X^ϵ_0   =   x_0,Y^ϵ_0=y_0


We assume non-degeneracy of the diffusion coefficients as well 𝒞^1 smoothness and boundedness of the drift and diffusion coefficients. Moreover, we assume that δ=δ(ϵ)↓0 such that ϵ/δ↑∞ as ϵ↓0.  (W_t, B_t) is a 2κ-dimensional standard Wiener process. We assume that for each fixed x∈ℝ^m,  b(·,γ), c(x,·,γ),σ(x,·,γ),f(·,γ), g(x,·,γ), τ_1(·,γ) and τ_2(·,γ) are stationary and ergodic random fields in an appropriate probability space (Γ,𝒢,ν) with γ∈Γ.


Notice that if we choose b(y,γ)=f(y,γ)=-∇_yQ(y,γ)  for a periodic function Q(·), c(x,y,γ)=-∇_xV(x),  σ(x,y,γ)=τ_1(y,γ)=√(2D) and τ_2(y,γ)=0, and set y_0=x_0/δ, we then get the Langevin equation (<ref>). In particular, if we make these choices, then we simply have   Y^ϵ_t=X^ϵ_t/δ and the model can be interpreted as  diffusion in the rough potential ϵ Q(x/δ,γ)+ V(x), where the roughness is dictated by  Q. In general, Q may not be modelled as a periodic function. One may model Q as a random field; see the simulation study in Subsection <ref>.




 §.§ Description of the random environment

The large deviations and importance sampling results for (<ref>), see <cit.>, are true under certain assumptions on the random medium that we recall here for convenience.  We assume that there is a group of measure preserving transformations {τ_y, y∈ℝ^d-m} acting ergodically on Γ that is defined as follows.




  * τ_y preserves the measure, namely ∀ y∈ℝ^d-m
and ∀ A∈𝒢 we have ν(τ_yA)=ν(A).


  * The action of {τ_y: y∈ℝ^d-m} is ergodic, that is if
A=τ_yA for every y∈ℝ^d then ν(A)=0 or 1.


  *  For every measurable function f on (  Γ, 𝒢,
ν), the function (y,γ)↦ f(τ_yγ) is
measurable on (  ℝ^d-m×Γ, 𝔹(ℝ
^d-m)⊗𝒢).



Let ϕ̃ be a square integrable function in Γ and define the operator T_yϕ̃(γ)=ϕ̃(τ_yγ). The operator T_y·  is a strongly continuous group of unitary maps in
L^2(Γ), see <cit.>. Denote by D_i the infinitesimal generator  of T_y in the direction
i, which is a closed and densely defined generator, see <cit.>.

In order to guarantee that the involved functions are ergodic and stationary random fields on ℝ^d-m,
for ϕ̃∈ L^2(Γ), let us define the operator ϕ(y,γ)=ϕ̃(τ_yγ).  Similarly, for a measurable function ϕ̃:ℝ^m×Γ↦ℝ^m
we consider the (locally) stationary random field (x,y) ↦ϕ̃(x,τ_yγ)=ϕ(x,y,γ). Then, it is guaranteed that ϕ(y,γ) (respectively ϕ(x,y,γ)) is a stationary (respectively locally stationary) ergodic random field.


The coefficients, b,c,σ,f,g,τ_1,τ_2 of (<ref>) are defined through this procedure and therefore are guaranteed to be ergodic and stationary random fields. For example in the case of the c drift term,  we start with an L^2(Γ) function c̃(x,γ) and we define the corresponding coefficients  via the relation c(x,y,γ)=c̃(x,τ_yγ).

For every γ∈Γ, let us the operator

    ℒ^γ=f(y,γ)∇_y·+tr[(
    τ_1(y,γ)τ^T_1(y,γ)+τ_2(y,γ)τ^T_2(y,γ))∇_y∇_y·]

which is the infinitesimal generator of a Markov process, say Y_t,γ. Using the Markov process Y_t,γ, we can define the so-called environment process, see <cit.>, denoted by γ_t. The environment process γ_t has continuous
transition probability densities with respect to the d-dimensional Lebesgue
measure, see <cit.>, and is defined by the equations

    γ_t     =τ_Y_t,γγ
    γ_0     =τ_y_0γ



The infinitesimal generator of the Markov process γ_t is given by

    L̃=f̃(γ)D·+tr[ ( τ̃_1(γ)τ̃_1^T(γ)+τ̃_2(γ)τ̃_2^T(γ))D^2·].


In order to simplify the presentation, let us assume that the operator  L̃ is in divergence form. In particular, let us set f̃(γ)=-D
Q̃(γ) and τ̃_1(γ)=√(2D)θ=constant and τ̃_2(γ)=√(2D)√(1-θ^2)=constant.

Then, we can write the unique ergodic invariant measure for the environment process {γ_t}_t≥ 0 in closed form; see <cit.> for more general case which is not necessarily restricted to the gradient case. Denote by 𝔼^ν  the expectation
operator with respect to the measure ν. Then , the  measure π(dγ) defined  on (Γ,𝒢) by

    π(dγ)≐m̃(γ)/𝔼^νm̃(·
    )ν(dγ),  with m̃(γ)=exp[-Q̃(γ)/D].

is the unique ergodic invariant measure for the environment process {γ_t}_t≥0.

Next, we need to define the equivalent to the cell problem in the case of periodic coefficients, also known as the macroscopic problem in the homogenization theory. To do so, we first define ℋ^1=ℋ^1(ν) to be the Hilbert space  equipped with the
inner product

    (f̃,g̃)_1=∑_i=1^d(D_if̃,D_ig̃).


 Let us consider ρ>0 and consider the following problem on Γ

    ρχ̃_ρ-L̃χ̃_ρ=b̃.


Under the condition b̃∈ L^2(ν) with ‖b̃‖_ℋ^-1<∞,  Lax-Milgram lemma, see <cit.>, guarantees that
equation (<ref>) has a unique weak solution in the abstract Sobolev space
ℋ^1 or equivalently in ℋ^1(π). At this point, we note that in the periodic case one also considers (<ref>), but one can then take ρ=0 given that b averages to zero when is integrated against the invariant measure π. However, in the random case, (<ref>) with ρ=0 does not necessarily have a well defined solution (even if b averages to zero when is integrated against the invariant measure π), see for example <cit.>.

In the general random case, we consider the equation with ρ>0 and in the end, the homogenization theorem is proven by taking appropriate sequences ρ=ρ(ϵ) such that ρ(ϵ)↓ 0 as ϵ↓ 0. Taking ρ↓ 0 is allowed by the following well known properties of the solution to (<ref>), (see <cit.>),


  * There is a constant K that is independent of ρ such that

    ρ𝔼^π[  χ̃_ρ(·)]  ^2
    +𝔼^π[  Dχ̃_ρ(·)]  ^2≤ K


  * χ̃_ρ has an ℋ^1
strong limit, i.e., there exists a χ̃_0∈ℋ^1(π)
such that

    lim_ρ↓0‖χ̃_ρ(·)-χ̃
    _0(·)‖ _1=0   and   lim_ρ↓0ρ𝔼^π[  χ̃_ρ
    (·)]  ^2=0.






 §.§ Large deviations and importance sampling theory for diffusion in random environments.

Now that we have defined the random environment and explained its properties, let us review the related large deviations and importance sampling theory from <cit.>. Set for notational convenience ξ̃=Dχ̃_0.



 Let {(X^ϵ,γ,Y^ϵ,γ),ϵ>0} be, for fixed γ∈Γ, the unique strong
solution to (<ref>). Assume non-degeneracy of the diffusion coefficients as well as 𝒞^1 smoothness and boundedness of the drift and diffusion coefficients. Consider the regime where ϵ,δ↓ 0 such that ϵ/δ↑∞. Then, {X^ϵ,γ,ϵ>0} converges in probability, almost surely with respect to the random environment  γ∈Γ, as ϵ,δ↓ 0 to the solution of the ODE

    dX̅_t=r(X̅_t)dt

and satisfies, almost surely with respect to γ∈Γ, the large deviations principle with rate function

    S_t_0T(ϕ)=
    1/2∫_t_0^T(ϕ̇_s-r(ϕ_s))^Tq^-1(ϕ_s
    )(ϕ̇_s-r(ϕ_s))ds    if ϕ∈𝒜𝒞
    ([t_0,T])   and ϕ_t_0=x_0
    
    +∞   otherwise.

where

    r(x)   =lim_ρ↓0𝔼^π[  c̃(x,·) +Dχ̃_ρ(·) g̃(x,·)]  =𝔼^π[c̃(x,·)+ξ̃
    (·)g̃(x,·)]
    
    q(x)   =lim_ρ↓0𝔼^π[  (σ̃(x,·)+Dχ̃_ρ
    (·)τ̃_1(·))(σ̃(x,·)+Dχ̃_ρ
    (·)τ̃_1(·))^T+(Dχ̃_ρ
    (·)τ̃_2(·))(Dχ̃_ρ
    (·)τ̃_2(·))^T]
       =𝔼^π[  (σ̃(x,·)+ξ̃(·)τ̃_1(·))(σ̃(x,·)+ξ̃
    (·)τ̃_1(·))^T+(ξ̃
    (·)τ̃_2(·))(ξ̃
    (·)τ̃_2(·))^T]

and ρ=ρ(ϵ)=δ^2/ϵ.


 Notice that the coefficients r(x) and q(x) are obtained by homogenizing (<ref>) by taking δ↓ 0 with ϵ fixed. The  form of the action functional can be recognized as the one that would come up when considering large deviations for  the homogenized system. This is also implied by the fact that δ goes to zero faster than ϵ, since  ϵ/δ↑∞.

We also remark here that if b=0, then χ_ρ=0. In this case  r(x),q(x) take the simplified forms r(x)=𝔼^π[c̃(x,·)] and q(x)=𝔼^π[  σ̃(x,·) σ̃(x,·)^T].

Turning now to importance sampling, given controls u_1 and u_2 one considers the controlled dynamics under the importance sampling measure ℙ̅

    dX̅^ϵ_s   =   [  ϵ/δb(Y̅^ϵ_s,γ)+c(  X̅^ϵ_s
    ,Y̅^ϵ_s,γ)+σ(  X̅_s^ϵ,Y̅_s^ϵ,γ)  u_1(s)]   dt+√(ϵ)
    σ(  X̅^ϵ_s,Y̅^ϵ_s,γ)
    dW̅_s, 
    
    dY̅^ϵ_s   =   1/δ[  ϵ/δf(Y̅^ϵ_s,γ)  +g(  X̅^ϵ_s
    ,Y̅^ϵ_s,γ)+τ_1(Y̅^ϵ_s,γ)u_1(s)+
    τ_2(Y̅^ϵ_s,γ)u_2(s)]   dt
          +√(ϵ)/δ[
    τ_1(Y̅^ϵ_s,γ)
    dW̅_s+τ_2(Y̅^ϵ_s,γ)dB̅_s],
    X̅^ϵ_t_0   =   x_0,Y̅^ϵ_t_0=y_0


where (v_1(s), v_2(s)) denote the first and second component of the control

    u(s,X̅^ϵ_s,Y̅^ϵ_s)=(u_1(s,X̅^ϵ_s,Y̅^ϵ_s), u_2(s,X̅^ϵ_s,Y̅^ϵ_s)).


Then, for a given cost function h(x), under ℙ̅

    Δ^ϵ,γ(t_0,x_0,y_0)=exp{  -1/ϵh(X̅^ϵ_T)}dℙ/dℙ̅(X̅^ϵ, Y̅^ϵ),

is an unbiased estimator for
𝔼[exp{  -1/ϵh(X^ϵ_T)}].

Consider next the Hamiltonian

    H(x,p)=<r(x), p>-1/2‖ q^1/2(x)p‖^2

 with r(x),q(x) the coefficients defined in Theorem <ref> and consider the HJB equation associated to this Hamiltonian, letting U̅(t,x) be a smooth subsolution to it (analogously to Section <ref> with r(x) and q(x) in place of -∇ V(x) and Γ respectively). Then, the following theorem guarantees at least logarithmic asymptotically good performance.


 Let {(X^ϵ_s, Y^ϵ_s),ϵ>0} be
the solution to (<ref>) for s∈[t_0,T] with initial point (x_0,y_0) at time t_0.  Consider a non-negative, bounded and continuous
function h:ℝ^m↦ℝ. Let U̅(s,x) be a subsolution to the associated HJB equation that has continuous derivatives up to
order 1 in t and order 2 in x, and the first and second
derivatives in x are uniformly bounded. Assume non-degeneracy of the diffusion coefficients as well 𝒞^1 smoothness and boundedness of the drift and diffusion coefficients. In the general case where b≠ 0, consider ρ>0 and define the (random) feedback control u_ρ(s,x,y,γ)=(u_1,ρ(s,x,y,γ), u_2,ρ(s,x,y,γ)) by

 
    u_ρ(s,x,y,γ)=(-(σ+Dχ_ρτ_1)^T(x,y,γ)∇_xU̅(s,x), -(Dχ_ρτ_2)^T(y,γ)∇_xU̅(s,x))


Then for ρ=ρ(ϵ)=δ^2/ϵ↓ 0 we have that almost surely in γ∈Γ

    lim inf_ϵ→0-ϵln Q^ϵ,γ(t_0,x_0,y_0;u_ρ(·))≥
    G(t_0,x_0)+U̅(t_0,x_0).


If b=0, then set u(s,x,y,γ)=(-σ^T(x,y,γ)∇_xU̅(s,x), 0) and (<ref>) holds with u_ρ(·)=u(·).














 §.§ A simulation study

Consider for instance the case of Example <ref>

    dX^ϵ,δ_t=-∇
    V^ϵ(X^ϵ,δ_t,X^ϵ,δ_t/δ)dt+√(2ϵ)dW_t,

where the potential function
V^ϵ(x,x/δ)=ϵ Q(x/δ)+V(x).
Q(y) is a stationary ergodic random field on a probability space (𝒳,𝒢,ν). We may consider for instance V(x)=1/2x^2 and

    Q(y) mean zero Gaussian with E^ν[Q(x)Q(y)]=exp[-|x-y|^2]




Making the connection with (<ref>), the fast Y motion essentially is Y=X/δ. Referring to Theorems <ref> and <ref> we have r(x)=-V^'(x)/(KK̂) and q=2/(KK̂) where K=E^ν[  e^-Q(z)], K̂=E^ν[  e^Q(z)]. Given a classical subsolution U̅, one expects that the corresponding change of simulation measure that guarantees at least asymptotic optimality, is based on the control
u̅(s,x,y,γ)=(-√(2)(1+∂χ(y,γ)/∂ y)U̅_x(s,x),0) where one can compute that the weight function is 1+∂χ(y,γ)/∂ y=e^Q(y,γ)/K̂. Note that in contrast to the periodic case, the control u is random in that it implicitly depends on γ∈Γ, via the random field Q(y,γ).


Assume that we want to estimate

    θ^ϵ,δ=P[X^ϵ,δ hits  1  before  0|X^ϵ,δ_0=0.1]


As in Subsection <ref>, we compare the asymptotical optimal change of measure with standard Monte Carlo, which corresponds to no change of measure, and with the importance sampling that corresponds to the change of measure based only on the homogenized problem, which ignores the macroscopic problem.
Based on 10^7 trajectories, we have the following simulation data


It is clear, that the importance sampling scheme based on the asymptotically optimal change of measure u̅(t,x,y,γ) outperforms the standard Monte Carlo estimator in which no change of measure is being done. It also outperforms,  the estimator based solely on the homogenized system, which ignores the local information characterized by solution to the macroscopic problem. Of course, this behavior is parallel to the behavior observed in the periodic case of Subsection <ref>. Additional simulation studies can be found in <cit.>.

In <cit.>, the interested reader can find further simulation studies in the case of the general model (<ref>) where one does not necessarily have the Y motion to be X/δ. However, we do point out that the theoretical results of <cit.> are valid for the system (<ref>) where the process (X^ϵ,Y^ϵ) has initial point (x_0,y_0) and both x_0 and y_0 are of order one as δ↓ 0. This is not exactly the same to the case where Y=X/δ, as then y_0=x_0/δ, which is no longer of order one as δ↓ 0. But, simulation studies, as the one presented in Table <ref>, indicate that the theoretical results should be also valid for the Y=X/δ case.



§ IMPORTANCE SAMPLING FOR METASTABLE MULTISCALE PROCESSES AND FURTHER CHALLENGES

In Section <ref> we elaborated on the effects of rest points and metastable dynamics on importance sampling schemes. The end conclusion was that extra care is needed when stable or unstable equilibrium points are in the domain of interest. In this case, asymptotic optimality is not enough in that asymptotically optimal schemes may not perform well in practice unless one goes to really small values of ϵ, in which case the events may be too rare to be of any practical interest. Then, in Section <ref> and <ref> we summarized the issues that come up in the design of asymptotically efficient importance sampling schemes when the dynamics have multiple scales.

In <cit.> we have systematically addressed the effects of rest points onto the design of importance sampling schemes and have identified what the main issues are. In <cit.>, we have suggested a potential provably appropriate remedy to the issue, by constructions as the ones mentioned in Section <ref>. The subsolution constructed there effectively yields a very good approximation to the zero variance change of measure. 
Even though, the constructions in <cit.> work provably well pre-asymptotically and asymptotically and do not degrade as parameters such as the time horizon T getting large, the performance in higher dimensions can be worse than the corresponding  performance in the lower-dimensional cases. While this is expected to be the case as the dimension gets larger, due to further approximations and simplifications that need to be made, there is a clear room for improvement here. This is part of ongoing work of the author and we refer the interested reader to <cit.> for some results in the infinitely dimensional small noise SPDE case.

Moreover, it is clear that the constructions of Sections <ref> and <ref> guarantee only asymptotic optimality. If in addition to multiscale dynamics one has to also face metastability, then, as it was seen in Section <ref>, theoretical asymptotic optimality is not sufficient for good numerical performance. One can of course combine the results of Section <ref> with those of Sections <ref> and <ref>. To be more precise, one can combine the results of <cit.> with those of <cit.>. In practice, one can just use the changes of measure as indicated in <cit.> that guarantee asymptotic optimality, but construct the subsolution U̅(t,x) as indicated in <cit.>. We plan to address this issue in more detail in a future work.













99


AsmussenGlynn2007
S. Asmussen and P. W. Glynn, Stochastic Simulation: Algorithms and Analysis, Springer, (2007).


BLP A. Bensoussan, J.L. Lions and G. Papanicolaou, Asymptotic
analysis for periodic structures, Vol 5, Studies in Mathematics and its
Applications, North-Holland Publishing Co., Amsterdam, (1978).

BanushkinaMeuwly2007
P. Banushkina and M. Meuwly, Diffusive dynamics on multidimensional rough free energy surfaces, J. Chem. Phys. Vol. 127, (2007), 135101.



BlanchetGlynn
J.H. Blanchet and P. Glynn. Efficient rare-event simulation for the maximum of
heavy-tailed random walks. Ann. Appl. Prob., Vol. 18, (2008), pp. 1351-1378.

BlanchetGlynnLeder2011
J.H. Blanchet, P. Glynn and K. Leder,
On Lyapunov inequalities and subsolutions for efficient importance sampling, accepted ACM TOMACS, Vol. 22, Issue 3, (2012), Article No. 13.

BlanchetGlynnLiu
J.H. Blanchet, P. Glynn and J.C. Liu, Fluid heuristics, Lyapunov bounds and
efficient importance sampling for a heavy-tailed G/G/1 queue, QUESTA, Vol.113, (2007), pp. 57-99.

BlanchetLiu2008
J.H. Blanchet and J.C. Liu, State-dependent Importance Sampling for Regularly Varying Random Walks, Advances in Applied Probability, Vol. 40, (2008), pp 1104-1128.

BoyleBroadieGlasserman
P. Boyle, M. Broadie and P. Glasserman, Monte Carlo methods for security
pricing, Journal of Economic Dynamics and Control, Vol. 21, (1997), pp. 1257-1321.


BovierEckhoffGayrandKlein2004
 A. Bovier, M. Eckhoff, V. Gayrard and M. Klein, Metastability in reversible diffusion processes 1. Sharp estimates for capacities and exit times, J. Eur. Math. Soc. 6, (2004), pp. 399-424.

BovierGayrandKlein2005
A. Bovier, V. Gayrard and M. Klein. Metastability in reversible diffusion processes. 2. Precise estimates for small eigenvalues, J. Eur. Math. Soc. 7, (2005), pp. 69-99.











Cameron2012
M. Cameron, Finding the quasipotential for nongradient SDEs, Physica D: Nonlinear Phenomena, Vol. 241, Issue 18, (2012), pp. 1532–1550.



DayDarden1985
M. Day, T. Darden, Some regularity results on the Ventcel-Freidlin quasi-potential function, Appl. Math Opt. Vol. 13 (1985), pp. 259-282.


DemboZeitouni
A. Dembo and O. Zeitouni,
Large deviations techniques and applications, volume 38 of Applications of Mathematics,
Springer-Verlag, New York, 2nd ed., (1998).



DupuisLederWang2007
P. Dupuis, K. Leder and H. Wang, Large deviations and importance sampling
for a tandem network with slow-down, QUESTA, 57, (2007), pp. 71-83.

DupuisLederWang2009
P. Dupuis, K. Leder and H. Wang, Importance sampling for weighted serve-
the-longest-queue, Math. of Operations Research, Vol. 34 No. 3 (2009), pp. 642-660 .

DupuisSezerWang2007
P. Dupuis, A. Sezer and H. Wang, Dynamic importance sampling for queueing
networks, Annals of Applied Probability,  (2007), pp. 1306-1346.

DupuisSpiliopoulos
P. Dupuis and K. Spiliopoulos.
Large deviations for multiscale problems via weak convergence methods,
Stochastic Processes and their Applications, 122, (2012), pp. 1947-1987.


DupuisSpiliopoulosWang
P. Dupuis, K. Spiliopoulos and H. Wang.
Importance sampling for multiscale diffusions,
Multiscale Modeling and Simulation, Vol. 12, No. 1, pp. 1-27.

DupuisSpiliopoulosWang2
P. Dupuis, K. Spiliopoulos and H. Wang.
Rare Event Simulation in Rough Energy Landscapes.
2011 Winter Simulation Conference, pp. 504-515.


DupuisSpiliopoulosZhou2013
P. Dupuis, K. Spiliopoulos and X. Zhou, Escape from an equilibrium: importance sampling and rest points I,
 Annals of Applied Probability, Vol. 25, No. 5, (2015), pp. 2909–2958.

DupuisSpiliopoulos2014a
P. Dupuis and K. Spiliopoulos, Rare event simulation in the neighborhood of a rest point,
2014 Winter Simulation Conference, (IEEE, 2014), pp. 564-573.


DupuisWang
P. Dupuis and H. Wang, Importance sampling, large deviations and
differential games, Stochastics and Stochastics Reports,  76 (2004),
pp. 481-508.

DupuisWang2
P. Dupuis and H. Wang, Subsolutions of an Isaacs equation and efficient schemes of importance sampling, Mathematics of Operations Research, Vol.  32 (2007), pp. 723-757.




FlemingJames1992
W. H. Fleming and M. R. James, Asymptotic Series and Exit Time Probabilities,
Annals of Probability,  Vol. 20, No. 3, (1992), pp. 1369-1384.


FleimingSoner2006
W. H. Fleming and H. M. Soner,  Controlled Markov Processes and Viscosity Solutions, Springer, 2nd edition, 2006.

FreidlinWentzell88
M.I. Freidlin and A. D. Wentzell,
Random perturbations of dynamical systems, 2nd Edition,
Springer-Verlag, New York, (1988).


FS M. Freidlin and R, Sowers, A comparison of homogenization and
large deviations with applications to wavefront propagation,
Stochastic Process and Their Applications, Vol. 82, (1999), pp.
23–52.


Glasserman2004
P. Glasserman, Monte Carlo Methods in Financial Engineering, Springer-Verlag,
New York, 2004.

GlassermanKou1995
P. Glasserman and S. Kou, Analysis of an important sampling estimator for tandem queues, ACM Trans. Modeling Comp. Simulation, Vol. 4,  (1995), pp. 22-42.


GlassermanWang1997
P. Glasserman and Y. Wang, Counter examples in importance sampling for large deviations probabilities, Annals of Applied Probability, Vol. 7, (1997), pp. 731-746.


GlynnIglehart1988
P.W. Glynn and D.L. Iglehart,  Simulation methods for queues: an overview.
Queueing Systems: Theory and Applications, Vol. 3, (1988), pp. 221-256.

GriffithsTavare1994
R.C. Griffiths and S. Tavaré, Simulating probability distributions in the coalescent, Theoretical Population Biology, Vol. 46, (1994), pp. 131-159.

Janke
W. Janke, Rugged Free-Energy Landscapes, Lecture Notes in Physics, Volume 736/2008, Springer, (2008).



KomorowskiLandimOlla2012T. Komorowski, C. Landim and S. Olla,
Fluctuations in Markov Processes: Time Symmetry and Martingale
Approximation, Springer, 2012.


KosyginaRezakhanlouVaradhan
E. Kosygina, F. Rezakhanlou and S. R. S. Varadhan, Stochastic Homogenization of Hamilton-Jacobi-Bellman Equations, Communications on Pure and Applied Mathematics, Vol. LIX,  (2006), pp.0001–0033.




Levine
R.D. Levine, Monte Carlo, maximum entropy and importance sampling, Chemical Physics,  228, (1998), pp. 255-264.



LifsonJackson S. Lifson and J.L. Jackson, On the self-diffusion of
ions in a polyelectrolyte solution, Journal of Chemical
Physics, Vol. 36, (1962), pp. 2410-2414.



MaierStein1993
R.S. Maier and D.L. Stein, Escape problem for irreversible systems, Physical Review E, Vol. 48, No. 2, (1993), pp. 931-938.

MaierStein1997
R.S. Maier and D.L. Stein, Limiting exit location distributions in the stochastic exit problem, SIAM Journal on Applied Mathematics,
Vol. 57 Issue 3, (1997), pp. 752–790.


MazonkaJarzynskiBlocki1998
O. Mazonka, C. Jarzynski and J. Blocki, Computing probabilities of very rare
events for Langevin processes: A new method based on importance sampling.
Nuclear Physics A, 641, (1998), pp. 335-354.



Olla1994S. Olla, Homogenization of Diffusion Processes in Random
Fields, 1994, available at www.ceremade.dauphine.fr/∼ olla/lho.ps.






Osada1983H. Osada, Homogenization of diffusion processes with
random stationary coefficients,  Probability Theory and Mathematical
Statistics, Lecture Notes in Mathematics, Springer, Vol. 1021, (1983), pp 507-517.


PS G.A. Pavliotis and A.M. Stuart, Multiscale methods: Averaging and
Homogenization, Springer, 2007.

PapanicolaouVaradhan1982G. Papanicolaou and S.R.S. Varadhan,
Boundary value problems with rapidly oscillating random coefficients,
Colloquia Mathematica Societatis Janos Bolyai 27, Random Fields,
Esztergom (Hungary) 1979, North Holland, (1982), pp. 835-873.


SalinsSpiliopoulos2016
M. Salins and K. Spiliopoulos, Rare event simulation via importance sampling for linear SPDE's, Stochastics and Partial Differential Equations: Analysis and Computations, to appear, (2017).


Siegmund1976
D. Siegmund, Importance sampling in the Monte Carlo study of sequential tests, Annals of Statistics, 4, (1976), pp. 673-684.


Spiliopoulos2013a
K. Spiliopoulos,
Large deviations and importance sampling for systems of slow-fast motion,
  Applied Mathematics and Optimization,
Vol 67,  (2013), pp. 123-161.



Spiliopoulos2013b
K. Spiliopoulos,
 Non-asymptotic performance analysis of importance sampling schemes for small noise diffusions,
Journal of Applied Probability,  Vol. 52, (2015), pp. 1-14.

Spiliopoulos2013c
K. Spiliopoulos,
 Quenched Large Deviations for Multiscale Diffusion Processes in Random Environments,
 Electronic Journal of Probability,Vol. 20, (2015), no. 15, pp. 1�29.

Spiliopoulos2015
K. Spiliopoulos,
 Rare event simulation for multiscale diffusions in random environments,
 SIAM Multiscale Modeling and Simulation,  (2015), Vol. 13, No. 4, pp. 1290–1311.



VandenEijndenWeare
E. Vanden-Eijnden and J. Weare, Rare event simulation with vanishing error for small noise diffusions, Communications in Pure and Applied Mathematics,  Vol. 65(12), (2012), pp. 1770-1803.

Varadhan
S. R. S. Varadhan,
Large deviations and applications, volume 46 of CBMS-NSF
  Regional Conference Series in Applied Mathematics,
Society for Industrial and Applied Mathematics (SIAM), Philadelphia,
  PA, (1984).

VielPatelNiyazWhaley2002
A. Viel, M.V. Patel, P. Niyaz and K. Whaley, Importance sampling in rigid
body diffusion Monte Carlo, Computer Physics Communications, Vol. 145,
(2002), pp. 24-47.







ZuckermanWoolf
D. Zuckerman and T. Woolf, Efficient dynamic importance sampling of rare events in one dimension, Physical Review E., 63:016702, (2000), pp. 1-10.



Zwanzig R. Zwanzig, Diffusion in a rough potential, Proc. Natl.
Acad. Sci. USA, Vol. 85, (1988), pp. 2029-2030.




