Importance Sampling of Rare Events in Chaotic Systems

TL;DR

This paper introduces a novel Monte Carlo sampling method for efficiently sampling rare trajectories in chaotic dynamical systems, enabling the computation of distributions of rare events like Lyapunov exponents and escape times.

Contribution

It develops a Metropolis-Hastings Monte Carlo framework tailored for chaotic systems to sample exponentially rare events with polynomial computational effort.

Findings

Successfully computes distributions of finite-time Lyapunov exponents.

Efficiently samples escape times in transient-chaos scenarios.

Provides open-source software for implementation.

Abstract

Finding and sampling rare trajectories in dynamical systems is a difficult computational task underlying numerous problems and applications. In this paper we show how to construct Metropolis- Hastings Monte Carlo methods that can efficiently sample rare trajectories in the (extremely rough) phase space of chaotic systems. As examples of our general framework we compute the distribution of finite-time Lyapunov exponents (in different chaotic maps) and the distribution of escape times (in transient-chaos problems). Our methods sample exponentially rare states in polynomial number of samples (in both low- and high-dimensional systems). An open-source software that implements our algorithms and reproduces our results can be found in https://github.com/jorgecarleitao/chaospp

Equations

$${\mathbf{x}}_{t+1}=\mathbf{F}({\mathbf{x}}_t)={\mathbf{F}}^{t+1}({\mathbf{x}}_0),$$

$${\mathbf{x}}_t^{\prime }-{\mathbf{x}}_t={\mathbf{F}}^t({\mathbf{x}}^{\prime })-{\mathbf{F}}^t(\mathbf{x}).$$

$${\mathbf{x}}_t^{\prime }-{\mathbf{x}}_t=J_t(\mathbf{x})\cdot \mathbf{h}+\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n\frac{{\partial }^2{\mathbf{F}}^t(\mathbf{x})}{\partial x_i\partial x_j}{h}_i{h}_j+O(|\mathbf{h}{|}^3)$$

$$D(\mathbf{x},\mathbf{h},t)\equiv J_t(\mathbf{x})\cdot \mathbf{h}=\left(\prod _{i=t-1}^{0}J\left({\mathbf{x}}_i\right)\right)\cdot \mathbf{h}={\mathbf{h}}_t$$

$${\mathbf{h}}_0=\mathbf{h};{\mathbf{h}}_{i+1}=J\left({\mathbf{x}}_i\right)\cdot {\mathbf{h}}_i.$$

$${\lambda }_t(\mathbf{x})=\frac{\log ({\mu }_1)}{t},$$

$${\mathbf{x}}_t^{\prime }-{\mathbf{x}}_t={\delta }_0{e}^{{\lambda }_t(\mathbf{x})t}.$$

$${\lambda }_t(\mathbf{x})=\frac{1}{t}\sum _{i=0}^{t-1}\log |\frac{dF({\mathbf{x}}_i)}{dx}|.$$

$$P(E)=\int \delta (E-t_o{\lambda }_{t_o}(\mathbf{x}))U(\mathbf{x})\mathbf{d}\mathbf{x}.$$

$$W({\lambda }_{t_o})\equiv {\int }_{\Gamma }W(\mathbf{x})\delta ({\lambda }_{t_o}-{\lambda }_{t_o}(\mathbf{x}))U(\mathbf{x})\mathbf{d}\mathbf{x}$$

$$\mathbb{E}\left[W_f|{\lambda }_{t_o}\right]=\frac{1}{P({\lambda }_{t_o})}{\int }_{\Gamma }{W}_f(\mathbf{x})\delta ({\lambda }_{t_o}-{\lambda }_{t_o}(\mathbf{x}))U(\mathbf{x})\mathbf{d}\mathbf{x},$$

$$t_e(\mathbf{x})\equiv \min \{t:{\mathbf{F}}^t(\mathbf{x})\in \Lambda \}.$$

$$P(t_e)={\int }_{\Gamma }{\delta }_{t_e,t_e(\mathbf{x})}U(\mathbf{x})\mathbf{d}\mathbf{x},$$

$$W(t_e)\equiv {\int }_{\Gamma }W(\mathbf{x}){\delta }_{t_e,t_e(\mathbf{x})}U(\mathbf{x})\mathbf{d}\mathbf{x}.$$

$$\mathbb{E}\left[W|t_e\right]\equiv \frac{1}{P(t_e)}{\int }_{\Gamma }W(\mathbf{x}){\delta }_{t_e,t_e(\mathbf{x})}U(\mathbf{x})\mathbf{d}\mathbf{x}.$$

$$W(E)={\int }_{\Gamma }W(\mathbf{x})\delta (E,E(\mathbf{x}))U(\mathbf{x})\mathbf{d}\mathbf{x}.$$

$$\overline{W(E)}\equiv \frac{1}{m}\sum _{i=1}^m\delta (E,E(\mathbf{x}))W({\mathbf{x}}_i)$$

$$ϵ(E)\propto \frac{1}{\sqrt{mG(E)}}.$$

$$\overline{W(E)}\equiv \frac{1}{m}\sum _{i=1}^m\delta (E,E({\mathbf{x}}_i))W({\mathbf{x}}_i)\frac{U({\mathbf{x}}_i)}{\pi ({\mathbf{x}}_i)},$$

$$ϵ(E)\propto \frac{1}{\sqrt{mG(E)\pi (E)}}.$$

$$m(E)=mG(E)\pi (E).$$

$$ϵ(E)\propto \sqrt{\frac{1+2T(E)}{mG(E)\pi (E)}},$$

$$P({\mathbf{x}}^{\prime }|\mathbf{x})=g({\mathbf{x}}^{\prime }|\mathbf{x})a({\mathbf{x}}^{\prime }|\mathbf{x}),$$

$$a({\mathbf{x}}^{\prime }|\mathbf{x})=\min \left(1,\frac{g(\mathbf{x}|{\mathbf{x}}^{\prime })}{g({\mathbf{x}}^{\prime }|\mathbf{x})}\frac{\pi ({\mathbf{x}}^{\prime })}{\pi (\mathbf{x})}\right).$$

$$\pi (\mathbf{x})=\pi (E(\mathbf{x}))\propto e^{-\beta E(\mathbf{x})}.$$

$$m(E)\propto mG(E)e^{-\beta E}\propto m{e}^{-\beta E+S(E)}.$$

$$\pi (\mathbf{x})=\pi (E(\mathbf{x}))\propto \frac{1}{G(E(\mathbf{x}))},E\in [E_{\min },E_{\max }],$$

$$m(E)=\text{const.}$$

$$\tau (N)\sim N^2.$$

$$\tau (t_e^{\max })\sim (t_e^{\max }{)}^2.$$

Full text

1. I Introduction

2. II Review of problems

3. II.1 Variability of trajectories’ chaoticity

4. II.2 Transient chaos

5. II.3 Summary: numerical problems in the study of rare events in chaotic systems

6. III Review of methods

7. III.1 Importance Sampling

8. III.2 Metropolis-Hastings algorithm

9. III.3 Sampling distribution

10. III.4 Characterisation of the efficiency

11. III.5 Summary: challenges for the application of Monte Carlo methods to chaotic systems

12. IV Proposal distribution

13. IV.1 Aim of the proposal distribution

14. IV.2 Propose correlated trajectories

15. IV.3 Guarantee local proposals

16. IV.4 Summary: how to propose

17. V Numerical Tests

18. V.1 Finite-time Lyapunov exponent

19. V.2 Transient chaos

20. VI Conclusions

21. VI.1 Summary of results

22. VI.2 Comparison to previous results

23. VI.3 Discussion

24. VII Appendices

25. VII.1 Appendix 1: Dynamical Systems

26. VII.2 Appendix 2: Efficiency of the uniform proposal

27. VII.3 Appendix 3: Simplified proposals

28. VII.4 Appendix 4: Algorithmic description

I Introduction

Extreme events play a crucial role in our society. Landslides, floods, meteorite collisions, solar flares, earthquakes are all events that are rare but often lead to catastrophic consequences to our well being Albeverio et al. (2006). Science often studies extreme events by recreating the process that generates them sufficiently many times. While in some cases the process can be reproduced experimentally, in many others (e.g., astronomical and weather events) the best we can do is to simulate physical models.

Physical models often contain non-linearities in the equations of motion that amplify fluctuations and give origin to extreme events. A famous effect of non-linearity is chaos, the extremely sensitive dependence on initial conditions. Chaos has a fascinating history dating back to the seminal work of Poincaré in 1890 Holmes (1990) and is today a well established field of research with applications in Biology, Geology, Economy, Chemistry, and Physics Motter and Campbell (2013); Ott (1993). Chaotic dynamics often hinders our ability to study the evolution of the system analytically, e.g. it forbids describing trajectories in a closed formula. In this situation, again, often the best we can do is to numerically simulate the model in a computer, a paradigm popularized by Lorenz since the 1960s Lorenz (1963); Motter and Campbell (2013). Extreme events are then studied statistically, over an ensemble of initial conditions.

In this paper we introduce a framework for performing numerical simulations in chaotic systems in such a way that rare trajectories are generated more likely than if they would be generated by chance. This is an importance-sampling Newman and Barkema (2002); Robert and Casella (2005) (or rare-event-simulation Bucklew (2004)) strategy that builds on methods that have proven successful in the characterization of rare configurations in various problems. The distinguishing feature of our framework is that it considers general classes of observables Lucarini et al. (2016) in deterministic chaotic Ott (1993); Lai and Tél (2010) systems, being therefore able to find and sample initial conditions leading to extreme events in different problems. We apply our framework to two traditional problems in the numerical exploration of chaotic systems:

• •

trajectories with high or low finite-time Lyapunov exponents;

• •

long-living trajectories in open systems.

The results shown in this paper are general, but simulations are done in simple dynamical systems (time-discrete maps with up to 16 dimensions). Our goal is to illustrate the generic computational challenges of sampling rare events in chaotic systems, in line with the tradition that simple systems often possess the basic mechanisms responsible for extreme events.

The importance of rare trajectories in chaotic systems, including methods designed to find and sample them, have been studied through different perspectives Sweet et al. (2001); Bollt (2005); Lai and Tél (2010); Dellago et al. (2002); Bolhuis et al. (2002); Tailleur and Kurchan (2007); Kitajima and Iba (2011); Laffargue et al. (2013); Patra and Keshavamurthy (2015); Wouters and Bouchet (2016); Lucarini et al. (2016). For example, the method stagger and dagger, used to find long-living trajectories of transiently chaotic systems, has been an important tool to characterize chaotic saddles Sweet et al. (2001); Lai and Tél (2010). Lyapunov weighted dynamics is a population Monte Carlo method that was successfully applied to find trajectories with atypically low or high chaoticity Tailleur and Kurchan (2007); Geiger and Dellago (2010). Similar genealogical particle analysis was applied to compute rare events in simple climate models Wouters and Bouchet (2016). Another powerful and widely used method within trajectory sampling is transition path sampling Dellago et al. (2002); Bolhuis et al. (2002), that has been used to sample rare trajectories (e.g. chemical reactions), typically influenced by thermal noise Bolhuis et al. (2002). These are typically trajectories that transit from one stable configuration to another stable configuration Bolhuis et al. (2002); Grünwald et al. (2008). These different methods achieve their goal through different, often ingenious, solutions. Here we aim to construct methods that are applicable to different classes of problems and to quantitatively understand the impact of different parameters and choices on the efficiency of the algorithm. We then show that some of the existing methods can be derived from our construction under suitable approximations.

Our framework relies on Metropolis-Hastings (MH) importance sampling Newman and Barkema (2002); Robert and Casella (2005), a well established numerical technique that has been used in statistical physics to study rare events since the 1950s. MH produces a random walk $\boldsymbol{x}\rightarrow\boldsymbol{x}^{\prime}$ in the phase-space of the system that generates initial conditions leading to extreme events (rare states) more often than they would be found by chance, consequently reducing the computational cost associated with their rareness. The flexibility of MH is confirmed by its success in numerous fields in Physics, Chemistry, Finance, among many others Newman and Barkema (2002); Robert and Casella (2005). While transition path sampling Bolhuis et al. (2002), Lyapunov weighted dynamics Tailleur and Kurchan (2007), and genealogical methods Wouters and Bouchet (2016) already use importance sampling in deterministic chaotic systems, three fundamental questions remains largely open: 1. can MH be used to systematically sample rare trajectories of (deterministic) chaotic systems? If yes, 2. how and 3. at what (computational) cost?

To answer these questions, we develop a systematic approach to apply MH in chaotic systems. The crucial step is the choice of the proposal distribution of the MH algorithm, the conditional probability of “trying” a state $\boldsymbol{x}^{\prime}$ in the phase-space of the system, given the current state $\boldsymbol{x}$. The question we have to answer in order for MH to work is: what proposal distribution guarantees that an observable of the trajectory starting at $\boldsymbol{x}^{\prime}$, e.g. its Lyapunov exponent, is similar to the same observable of the trajectory starting at $\boldsymbol{x}$? Our main contribution is a methodology to answer this question for a broad class of chaotic systems and observables. More specifically, we show how incorporating properties of trajectories of chaotic systems in the proposal is a necessary condition to obtain an efficient MH algorithm. This methodology allows to construct efficient MH algorithms to sample rare events in different problems and classes of chaotic systems. We expect the ideas and formalism presented here to find applications in other problems of chaotic systems and in the study of extreme events more generally. Therefore, we expect our results to be useful both to those studying extreme events in non-linear systems and to those studying numerical techniques.

The paper is organized as follows: Sec. II reviews traditional numerical problems in chaotic systems and shows how they can be formulated as a problem of sampling rare trajectories; Sec. III introduces the MH algorithm as a method to perform importance-sampling simulations in chaotic systems, and shows that naive approaches do not lead to an efficient MH; Sec. IV shows how to incorporate general features of chaotic systems, such as exponential divergence or self-similarity of some of its properties, in the proposal distribution in order to achieve an efficient MH; Sec. V presents numerical tests of the general framework that confirm the applicability of Monte Carlo algorithms to sample rare events in different classes of chaotic systems; Sec. VI summarizes our results and discusses its implications.

II Review of problems

We consider dynamical systems whose states $\boldsymbol{x}$ in a phase space $\Omega$, $\boldsymbol{x}\in\Omega\subset\mathbb{R}^{D}$, evolve in time from an initial condition $\boldsymbol{x}=\boldsymbol{x}_{0}$ according to

[TABLE]

where $\boldsymbol{F}(\boldsymbol{x})\in\Omega$ and $\boldsymbol{F}^{t}$ is $\boldsymbol{F}$ composed $t$ times, $\boldsymbol{F}^{t}(\boldsymbol{x})\equiv\boldsymbol{F}(\boldsymbol{F}(...\boldsymbol{F}(\boldsymbol{x})...))$. Such discrete-time systems can be obtained from continuous-time systems through a Poincare surface of section, a stroboscopic projection, or by the time discretization of the numerical integration of differential equations. We are also interested in the dynamics in the tangent space, which quantifies the divergence between two initial conditions Ott (1993). Specifically, the distance of a state displaced from $\boldsymbol{x}$ by $\boldsymbol{h}$, $\boldsymbol{x}^{\prime}=\boldsymbol{x}+\boldsymbol{h}$, to the original state $\boldsymbol{x}$ evolves in time according to

[TABLE]

Expanding $\boldsymbol{F}^{t}(\boldsymbol{x}^{\prime})$ around $\boldsymbol{F}^{t}(\boldsymbol{x})$ allows $\boldsymbol{x}^{\prime}_{t}-\boldsymbol{x}_{t}$ to be written as

[TABLE]

where $J_{t}(\boldsymbol{x})\equiv{\bf d}\boldsymbol{F}^{t}(\boldsymbol{x})/{\bf d}\boldsymbol{x}$ is the Jacobian matrix of $\boldsymbol{F}^{t}$, and $\partial^{2}\boldsymbol{F}^{t}(\boldsymbol{x})/(\partial x_{i}\partial x_{j})$ is the $(i,j)$ entry of the Hessian matrix of $\boldsymbol{F}$. The first term of Eq. 3 can be expanded using the derivative of the composition and be written as

[TABLE]

where $J\equiv J_{1}$ and $\boldsymbol{h}_{t}$ evolves in the tangent space according to

[TABLE]

For small $|\boldsymbol{h}|$, the growth of $\boldsymbol{x}^{\prime}_{t}-\boldsymbol{x}_{t}$ is characterized by the eigenvalues and eigenvectors of $J^{t}$ Cvitanovic et al. (2016). The largest finite-time Lyapunov exponent (FTLE) of a point $\boldsymbol{x}$ can be defined111This is sometimes denoted as the stability exponent. See Chapter 4 and 6 of Ref. Cvitanovic et al. (2016). as

[TABLE]

where $\mu_{1}$ is the real part of the largest eigenvalue of $J^{t}$. Thus, when at least one direction is unstable ($\mu_{1}>1$), $\boldsymbol{x}_{t}^{\prime}-\boldsymbol{x}_{t}$ increases exponentially with time, and at most by

[TABLE]

When the system is one dimensional, the “Jacobian matrix” is a single number, the product in Eq. 4 is a product of numbers, and the only “eigenvalue” is the result of this product. Thus, in this case Eq. 6 can be written as

[TABLE]

We now describe two computational problems in chaotic dynamical systems.

II.1 Variability of trajectories’ chaoticity

For a chaotic system, $\lambda_{L}\equiv\lambda_{t\rightarrow\infty}>0$ in Eq. 6. The variation of the (maximum) finite-time Lyapunov exponents $\lambda_{t}$ (FTLE) across different initial conditions characterizes the variation of chaoticity of the system, yielding a distribution of the FTLE, $P(\lambda_{t_{o}})$, or, equivalently, the distribution of $E\equiv t_{o}\lambda_{t_{o}}$:

[TABLE]

where $U(\boldsymbol{x})$ is a chosen probability distribution (e.g. uniform in the phase-space). The FTLE and its distribution was introduced in the 80s Grassberger and Kantz (1985); Grassberger et al. (1988); Sepúlveda et al. (1989) and has been used to study turbulent flows Beigie et al. (1993), Hamiltonian dynamics Amitrano and Berry (1992), chimera states Olmi (2015), characterize dynamical trapping Grassberger and Kantz (1985); Zaslavsky (2002); Szezech et al. (2005); Artuso and Manchein (2009); Manchein et al. (2012), among others Abarbanel et al. (1991, 1992). The distribution of FTLE is related to the generalized dimensions Beck and Schögl (1995) and often follows a large deviation principle, where $t_{o}$ is the extensive parameter Beck and Schögl (1995); Laffargue et al. (2013). In strongly chaotic systems, the distribution of FTLE is Gaussian Beck and Schögl (1995), whereas for intermittent chaos and other weakly chaotic systems the distribution is typically non-Gaussian Prasad and Ramaswamy (1999).

Figure 1 shows the phase-space dependency and distribution of the FTLE in one chaotic system. It contains characteristic features observed in strongly chaotic systems: $\lambda_{t}$ is a quantity that depends sensitively on the state $\boldsymbol{x}$, and its distribution $P(\lambda_{t})$ decays exponentially to zero on both sides. In the limit $t_{o}\rightarrow\infty$, $P(\lambda_{t_{o}})\rightarrow\delta(\lambda_{t_{o}}-\lambda_{L})$.

The tails of the distribution $P(E)$ play a significant role in the characterization of chaotic systems, as, for example, the higher moments of the distribution are related to higher $q$s in the generalized dimensions $D_{q}$ of the attractor Beck and Schögl (1995). Furthermore, the regions of the phase-space with small (large) finite-time Lyapunov exponent are associated with slow (fast) decay of correlations Sepúlveda et al. (1989), and their characterization has been used to get insight on whether the system is ergodic or not Amitrano and Berry (1992). Moreover, trajectories characterized by a low or high finite-time Lyapunov exponent can play a significant role in the dynamics of interfaces in chaotic flows Beigie et al. (1993) and others Tailleur and Kurchan (2007).

A typical analysis of the FTLE is to measure how a quantity, $W(\boldsymbol{x})$, depends on the FTLE $\lambda_{t_{o}}$ Sepúlveda et al. (1989); Amitrano and Berry (1992); Beigie et al. (1993). This requires estimating an integral of the form

[TABLE]

where $W(\boldsymbol{x})$ is the pre-selected quantity, $\Gamma$ is a pre-selected sampling region (often $\Gamma=\Omega$), and $U(\boldsymbol{x})$ is the weight attributed to $\boldsymbol{x}$ (often uniform, $U(\boldsymbol{x})=1/|\Gamma|$). Let us look for two examples. First, consider the problem of estimating the distribution $P(\lambda_{t_{o}})$, e.g. Refs. Sepúlveda et al. (1989); Amitrano and Berry (1992); Prasad and Ramaswamy (1999); Tailleur and Kurchan (2007). The traditional technique is to sample $M$ states $\boldsymbol{x}$ according to $U(\boldsymbol{x})=1/|\Gamma|$ and estimate $P(\lambda_{t_{o}})$ using the estimator $M_{\lambda}/M$, where $M_{\lambda}$ is the number of samples with $\lambda_{t_{o}}(\boldsymbol{x})\in[\lambda_{t_{o}},\lambda_{t_{o}}+\Delta\lambda_{t_{o}}]$. Formally, this corresponds to $W(\lambda_{t_{o}})$ with $W(\boldsymbol{x})=1$. The second example is retrieved from Ref. Sepúlveda et al. (1989). There, in order to evaluate the contribution of the algebraic region of the phase-space to the power spectra, the authors decomposed it in two terms corresponding to the power spectra of trajectories with low FTLE and trajectories with high FTLE. This required estimating the power spectra from a set of trajectories conditioned to an interval of $\lambda$s on the tails of the distribution. Associating the power spectrum ($S(f)$ in the ref.) with $W(\boldsymbol{x})=W_{f}(\boldsymbol{x})$ where $f$ is the frequency, the power spectra represented in Fig. 4 of Ref. Sepúlveda et al. (1989) corresponds to integrals (for different frequencies $f$) given by

[TABLE]

which contains integrals of the form of Eq. 10.

Numerically estimating the integral in Eq. 10 is challenging for two reasons: first, $P(\lambda_{t_{o}})$ decays with $t_{o}$ (lower panel of Fig. 1): the distribution of FTLE often follows a large deviation principle, $P(\lambda_{t_{o}}(\boldsymbol{x}))\propto\exp\left(t_{o}s(\lambda_{t_{o}})\right)$, where $s$ is intensive in respect to $t_{o}$ and is often concave Beck and Schögl (1995). Consequently, the traditional methodology of sampling states uniformly to find or sample states with increasing $t_{o}$ requires an exponentially high number of initial conditions. Second, the dependency of $\lambda_{t_{o}}(\boldsymbol{x})$ on $\boldsymbol{x}$ has multiple local minima and maxima (upper panel of Fig. 1). Such rough (fractal Richter (2016)) landscapes are known to challenge numerical techniques (e.g., simulations get trapped in local minima or maxima) Hansmann and Okamoto (1993).

The problem of finding and sampling states with high or low-$\lambda$ has been addressed in the literature with numerical techniques that go beyond traditional uniform sampling Tailleur and Kurchan (2007); Kitajima and Iba (2011); Laffargue et al. (2013). Such techniques have been successfully applied to find Tailleur and Kurchan (2007) and sample Kitajima and Iba (2011); Laffargue et al. (2013) states with extremal $\lambda$s in different chaotic systems. Refs. Tailleur and Kurchan (2007); Laffargue et al. (2013) use a population Monte Carlo canonical ensemble where $\lambda_{t}$ plays the role of the energy $E$ to find or sample states with high or low FTLE. The method computes stochastic trajectories that, from the numerical tests performed, are indistinguishable from (deterministic) trajectories of the system. Ref. Kitajima and Iba (2011) proposes a flat-histogram simulation to find high or low chaotic states by developing an observable to quantify the chaoticity of the state.

II.2 Transient chaos

The best known examples of chaotic systems have a fractal attractor in the phase space Ott (1993), the Lorenz attractor being the most prominent example Lorenz (1963). However, chaotic dynamics can appear also when the fractal invariant set in the phase space is not purely attracting, e.g. it may have stable and unstable directions (a saddle). Trajectories close to this set perform chaotic motion for an arbitrarily long (yet finite) time. This phenomenon appears in a variety of physical systems and is known as transient chaos Tél and Lai (2008); Altmann et al. (2013).

Numerical investigations of transiently-chaotic systems are computationally difficult because most trajectories quickly escape the vicinity of the chaotic saddle (on which the chaotic dynamics is properly defined) Sweet et al. (2001); Bollt (2005); Sala et al. (2016). More precisely, the escape time of a state $\boldsymbol{x}$ in a pre-selected region $\boldsymbol{x}\in\Gamma\subset\Omega$ is defined as the first passage time of its trajectory to a pre-selected exit set $\Lambda\subset\Omega$:

[TABLE]

Almost all trajectories start at $\Gamma$ and eventually leave when $\boldsymbol{x}_{t}\in\Lambda$, but they do so at different times $t_{e}$. Such variability is quantified by the distribution of escape time: the probability that a random initial condition $\boldsymbol{x}$ leaves at a given time $t_{e}$,

[TABLE]

where $\delta_{t_{e},t_{e}(\boldsymbol{x})}$ is the Kronecker delta and $U(\boldsymbol{x})$ is the probability density assigned to each state in $\Gamma$, which is often constant, $U(\boldsymbol{x})=1/\int_{\Gamma}{\bf d}\boldsymbol{x}\equiv 1/|\Gamma|$. Figure 2 shows the typical features of $t_{e}(\boldsymbol{x})$: it strongly depends on the initial state $\boldsymbol{x}$ and its distribution $P(t_{e})$ decays exponentially to zero (i.e., it has a constant escape rate $\kappa$).

There are two main numerical techniques to study transiently chaotic systems. The first is to find one long living state $\boldsymbol{x}$ with $t_{e}(\boldsymbol{x})\gg 1$ and to compute an average over states of this trajectory Lai and Tél (2010); Sweet et al. (2001). For large $t_{e}(\boldsymbol{x})$, the trajectory between times $[t_{e}(\boldsymbol{x})/2-t_{s},t_{e}(\boldsymbol{x})/2+t_{s}]$ where $t_{s}\ll t_{e}(\boldsymbol{x})/2$ is close to a trajectory on the chaotic saddle. When the saddle is ergodic, an average over this long trajectory corresponds to an average over the natural invariant density and therefore an average over these states characterizes invariant properties of the system. The second technique, which we focus in this work, is to compute averages over an ensemble $U(\boldsymbol{x})$ of initial conditions $\boldsymbol{x}$ that leave the system at time $t_{e}(\boldsymbol{x})=t_{e}$ Lai and Tél (2010); Altmann et al. (2013). For small $t_{e}$ observations depend on the particular initial density $U(\boldsymbol{x})$ (a point of interest in itself Altmann et al. (2013)), while for large $t$ they characterize invariant properties of the system (like in the previous technique, the states $\boldsymbol{F}^{t/2}(\boldsymbol{x})$ with $t\rightarrow\infty$ are independent samples of the natural invariant density).

A typical analysis within sampling transiently chaotic systems is to measure how a quantity, $W(\boldsymbol{x})$, changes with increasing escape time $t_{e}$. Numerically, this can be written as

[TABLE]

Let us enumerate 3 examples of computational problems that can be interpreted as numerical estimations of an integral of the form of $W(t_{e})$ in Eq. 14.

Compute the escape time distribution:

Numerically, $P(t_{e})$ is often computed by drawing states from $U(\boldsymbol{x})$ (e.g. uniform density), and counting the relative number of states that exited at escape time $t_{e}$ Lai and Tél (2010). This corresponds to $W(\boldsymbol{x})=1$ in which case Eq. 14 reduces to Eq. 13.

Compute generalized dimensions of the chaotic saddle:

The generalized dimensions are an important property of the chaotic saddle and its calculation is often performed by box counting Dhamala et al. (2001); Lai and Tél (2010); Ott (1993). Essentially, the phase-space is divided in equal, non-overlapping, and space-filling boxes $i=1,...,B(\varepsilon)$ of linear size $\varepsilon$ (intervals in 1 dimension, squares in 2, etc.) and the generalized dimension $D_{q}$ of exponent $q$ is proportional to $\log\sum_{i}^{B}(\varepsilon)\mu_{i}^{q}$, where $\mu_{i}$ is the fraction of points of the saddle that belong to the box $i$ Dhamala et al. (2001). Numerically, $\mu_{i}$ is estimated by first obtaining a set of points $\boldsymbol{x}_{j}$ in the saddle and then counting how many are in the particular box $i$. Such an estimate can be written as the expectation of an indicator function that tells whether a state in the saddle, $\boldsymbol{F}^{t_{e}/2}(\boldsymbol{x})$ for $t_{e}(\boldsymbol{x})\gg 0$, is inside the box $i$, $W(\boldsymbol{x})=\delta_{\boldsymbol{F}^{t_{e}/2}(\boldsymbol{x})\in i}$. This expectation can be written as a conditional expectation of $W(\boldsymbol{x})$ over states that leave at time $t_{e}$,

[TABLE]

Computing this essentially requires computing integrals of the form of $W(t_{e})$ in Eq. 14.

Compute the distribution of FTLE on the chaotic saddle:

The distribution of FTLE $P(\lambda)$ is another important property of the chaotic saddle Dhamala et al. (2001); Lai and Tél (2010). The FTLE $\lambda=\lambda_{t}$ for a fixed $t$ of an open system is computed for trajectories on the chaotic saddle. Like in the previous problem, each of these trajectories can be obtained by generating a state $\boldsymbol{x}$ according to $U(\boldsymbol{x})$ (e.g. uniformly distributed in $\Gamma$) that has a large escape time, $t_{e}(\boldsymbol{x})\gg 1$, and compute $\lambda_{t}(\boldsymbol{F}^{t_{e}/2}(\boldsymbol{x}))$ using Eq. 6 for a fixed $t$. The distribution of FTLE is then computed from an ensemble of these high-escape-time states by constructing different bins of the histogram $I_{\lambda}=[\lambda,\lambda+\Delta\lambda]$, and numerically compute the relative number of states in each bin. Formally, this equates to compute the expected number of states with a given escape time $t_{e}$ whose FTLE is in a bin, and thus corresponds to computing a conditional expectation of the form of Eq. 15 with $W(\boldsymbol{x})=W_{\lambda}(\boldsymbol{x})=\delta_{\lambda(\boldsymbol{F}^{t_{e}/2}(\boldsymbol{x}))\in I_{\lambda_{t}}}$.

II.3 Summary: numerical problems in the study of rare events in chaotic systems

To provide an unified treatment of the numerical problems in transient chaos and in computing FTLE of closed systems we use a common notation whenever possible. It indicates also how the methods can be generalized to different problems. Firstly, there is a quantity that we call ”observable” and denote by $E(\boldsymbol{x})$:

• •

$E(\boldsymbol{x})=t_{o}\lambda_{t_{o}}(\boldsymbol{x})$ in FTLE of closed systems

• •

$E(\boldsymbol{x})=t_{e}(\boldsymbol{x})$ in strongly chaotic open systems

Secondly, we use $\delta(E,E^{\prime})$ to denote both the Kronecker delta and the Dirac delta ($\delta(E-E^{\prime})$), for discrete and continuous case respectively. This is needed because, even though in practice we will always consider $E$ to be a discrete function due to binning of the histograms, formally the observable $E$ is either discrete ($t_{e}$ for maps) or continuous ($\lambda_{t}$ or $t_{e}$ for flows). Thirdly, the observable $E$ is a function of the phase-space of the system but typically depends also on an external quantity $N$ that parameterizes how high $E$ can be:

• •

$t_{o}$ in FTLE of closed systems

• •

(a pre-selected) maximum escape time, $t_{\max}$, in strongly chaotic open systems.

This parameter is important to us because, the higher it is, the rarer a state with the maximum or minimum possible $E$ is. In this notation, the two general problems presented in the previous sections can be summarized as follows:

• •

the analysis is conditioned to a projection of $\boldsymbol{x}$ into a one-dimensional observable $E(\boldsymbol{x})$.

• •

the probability distribution of $E$, $P(E)$, decays exponentially with increasing $E$.

• •

the focus of the analysis is on rare states in respect to $P(E)$: states with $E$ in one of the tails of $P(E)$.

• •

the rareness increases with $N$.

In transient chaos problems we are interested in states with increasing $N=t_{\max}$ that correspond to the tail of $P(E)=P(t_{e})$, $P(t_{\max})$. In computations of the FTLE in closed systems we are interested in states with increasing $N=t_{o}$ and on the tails of $P(E)$. These problems share two distinct computational problems: find rare states Sweet et al. (2001); Bollt (2005); Lai and Tél (2010); Tailleur and Kurchan (2007) and sample rare states Kitajima and Iba (2011); Laffargue et al. (2013); Leitão et al. (2013, 2014), which can be formalized as follows:

• •

Find rare states: minimize or maximize $E(\boldsymbol{x})$ over a constraining region $\Gamma$, $\boldsymbol{x}\in\Gamma\subset\Omega$, for increasing $N$.

• •

Sample rare states: compute the integral of a function $W(\boldsymbol{x})$ conditioned to a particular value of $E$ and for increasing $N$, over an ensemble of states $\boldsymbol{x}$ distributed according to $U(\boldsymbol{x})$ in a constraining region of the phase-space $\boldsymbol{x}\in\Gamma\subset\Omega$:

[TABLE]

The numerical challenges of these two numerical problems are common: the states $\boldsymbol{x}$ are exponentially difficult to find with increasing $N$ (or $E$), the function $E(\boldsymbol{x})$ contains multiple local and global minima embedded in fractal-like structures, and a potential high dimensionality $D$ of the phase-space (see Figs. 1 and 2). The goal of this paper is to develop a systematic approach to tackle the two numerical problems and these three challenges.

III Review of methods

III.1 Importance Sampling

As described in the previous section, the traditional methodology to compute an integral of the form of Eq. 16 is to draw $m$ samples $\{\boldsymbol{x}_{i}\}$ distributed according to $U(\boldsymbol{x})$ from the relevant region $\Gamma$ and approximate the integral $W(E)$ by the estimator

[TABLE]

where $\delta(E,E(\boldsymbol{x}))=1$ when $E(\boldsymbol{x}_{i})\in[E,E+\Delta E[$ and zero otherwise (when E is discrete, $\Delta E=1$). The relative distance of the estimator $\overline{W(E)}$ to $\mathbb{E}\left[W(E)\right]$ is quantified by the ratio $\epsilon(E)\equiv\sigma\left[\overline{W(E)}\right]/\mathbb{E}\left[W(E)\right]$ where $\sigma\left[\overline{W(E)}\right]^{2}\equiv\mathbb{E}\left[\overline{W(E)}^{2}\right]-\mathbb{E}\left[\overline{W(E)}\right]^{2}$. For the estimator in Eq. 17, this is given by

[TABLE]

where $G(E)$ is the density of states: the number of states per bin with an observable $E$.222$G(E)=P(E)$ when each state is equally weighted, $U(\boldsymbol{x})=1/|\Gamma|$. Therefore, the number of samples $m_{*}$ required to achieve a given precision $\epsilon_{*}$ for a given $E_{*}$ is $m_{*}(E_{*})\propto 1/G(E_{*})$. The critical problem in sampling rare states is that, because $G(E)$ decays exponentially with $E$, $m_{*}(E)$ increases exponentially with $E$.

Importance sampling techniques aim to improve this scaling by drawing samples from a distribution $\pi(\boldsymbol{x})\neq U(\boldsymbol{x})$ on the phase-space (e.g. non-uniformly in the phase-space) Robert and Casella (2005). Specifically, consider $m$ independent samples $\{\boldsymbol{x}_{i}\}$ drawn from $\pi(\boldsymbol{x})$, and the function $\pi(\boldsymbol{x})$ to depend only on $E$, $\pi(\boldsymbol{x})=\pi(E(\boldsymbol{x}))=\pi(E)$. Because the samples are not drawn from $U(\boldsymbol{x})$, the estimator in Eq. 17 would be biased. Importance sampling uses an unbiased estimator for $W(E)$ given by Robert and Casella (2005)

[TABLE]

which reduces to Eq. 17 when $\pi(\boldsymbol{x})=U(\boldsymbol{x})$. The advantage of importance sampling is that the relative error of the estimator in Eq. 19 is given by

[TABLE]

This is because, when sampling from $\pi(\boldsymbol{x})=\pi(E(\boldsymbol{x}))$, the expected number of samples with a given $E$, $m(E)$, is equal to

[TABLE]

Equation 20 implies that the function $\pi(E)$ can be chosen to favour states $\boldsymbol{x}$ with observable $E$ on the tails of $G(E)$ and therefore improve the precision of the estimator on these tails.

The standard deviations in Eqs. 18,20 were obtained assuming that the $m$ samples were independent. In traditional methodologies such as uniform sampling, this is the case. However, in the algorithms discussed below, it is not. Therefore, it is necessary to modify Eq. 20 for the case where the samples $\{\boldsymbol{x}_{i}\}$ are drawn from $\pi(\boldsymbol{x})$ and are also correlated. This modification is given by

[TABLE]

where $T(E)$ is the autocorrelation time Robert and Casella (2005), which increases with the correlation of the samples.

Efficient importance-sampling Newman and Barkema (2002); Robert and Casella (2005); Bucklew (2004) techniques have to address the following three steps:

choose a suitable $\pi(\boldsymbol{x})$ 2. 2.

have a method to generate samples from $\pi(\boldsymbol{x})$ 3. 3.

minimize the autocorrelation time $T$

A defining point in our method is our choice for the Metropolis-Hasting algorithm to address point 2, differently from Refs. Tailleur and Kurchan (2007); Wouters and Bouchet (2016) which address similar problems through a different choice for point 2 (cloning trajectories). Below we first discuss point 2, then point 1, and finally point 3.

III.2 Metropolis-Hastings algorithm

The Metropolis-Hastings (MH) algorithm asymptotically generates states $\boldsymbol{x}$ according to $\pi(\boldsymbol{x})$ using a a Markovian, ergodic and detailed balance random walk in the sampling region $\Gamma$ Robert and Casella (2005). This random walk is initialized from a random state $\boldsymbol{x}\in\Gamma$ and evolves to a new state $\boldsymbol{x}^{\prime}\in\Gamma$ with a transition probability $P(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ chosen such that asymptotically the states $\boldsymbol{x}$ are drawn according to $\pi(\boldsymbol{x})$. In Metropolis-Hastings, $P(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ is written as

[TABLE]

where $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ is the (proposal) distribution used to generate new states and $a(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ is the (acceptance) distribution used to select them. The random walk fulfills detailed balance because the acceptance probability is chosen as Robert and Casella (2005)

[TABLE]

Metropolis-Hastings algorithm is not the only way to achieve this. Another popular and alternative method to sample from $\pi(\boldsymbol{x})$ is population Monte Carlo, which instead of a random walk, uses multiple stochastic trajectories that are cloned and destroyed. See e.g. ref. Laffargue et al. (2013) for an application to the problem of the FTLE described above.

Algorithmically, the MH algorithm is implemented as follows: choose a region $\Gamma$ and a random initial condition $\boldsymbol{x}\in\Gamma$. Evolve the random walk in time according to:

Propose a state $\boldsymbol{x}^{\prime}$ drawn from $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})$; 2. 2.

Compute $a(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ replacing $\boldsymbol{x}$ and $\boldsymbol{x}^{\prime}$ in Eq. 23; 3. 3.

Generate a random number $r$ in $[0,1]$. If $r<a(\boldsymbol{x}^{\prime}|\boldsymbol{x})$, make $\boldsymbol{x}^{\prime}$ to be the new $\boldsymbol{x}$; 4. 4.

Store $\boldsymbol{x}$ and go to 1.

The set of sub-steps 1-4 brings the random walk from its current state $\boldsymbol{x}$ to the next state and it is called a Markov step. After a transient number of steps where the algorithm converges to the asymptotic distribution, the stored states $\boldsymbol{x}$ are (correlated) samples drawn from $\pi(\boldsymbol{x})$, and can be directly used in Eq. 19.

III.3 Sampling distribution

Canonical ensemble

A sampling distribution $\pi(\boldsymbol{x})$ often used is the canonical distribution Newman and Barkema (2002)

[TABLE]

In the context of rare states, the canonical ensemble is useful because the number of sampled states $m(E)$ in Eq. 21 becomes

[TABLE]

In particular, the maximum of $m$ is at $E^{*}$ solution of $\beta=dS/dE(E^{*})$. Therefore, $\beta$ tunes which value of the observable $E$ is sampled the most. For example, the Lyapunov weighted dynamics in Ref. Laffargue et al. (2013) uses this distribution (Eq. 15 of the ref. with $\alpha$ replaced by $\beta$).

Flat-histogram ensemble

Another distribution often used in the literature of Metropolis-Hastings Torrie and Valleau (1977); Swendsen and Wang (1986); Marinari and Parisi (1992); Sugita and Okamoto (1999); Lee (1993); Berg and Neuhaus (1991); Wang and Landau (2001); Viana Lopes et al. (2006) is the flat-histogram (or multicanonical), given by Berg and Neuhaus (1991); Wang and Landau (2001)

[TABLE]

for a given choice of $E_{\min},E_{\max}$ that defines the region of interest on the observable $E$. This is known as flat-histogram because, replacing Eq. 26 in Eq. 21 leads to a constant average number of samples on each $E$,

[TABLE]

Consequently, the dependence of the variance in Eq. 22 is only due to the autocorrelation $T(E)$, which implies that the computational cost to draw a state on the tail of $G(E)$ is no longer driven by the exponential decrease of $G(E)$, but by the computational cost to draw uncorrelated samples from $\pi(\boldsymbol{x})$.

The main limitation of the flat-histogram is that it requires knowing $G(E)$ in advance, which is very often unknown. The most well known, that we use here, is the Wang-Landau algorithm Wang and Landau (2001), that modifies the Metropolis-Hastings algorithm to a non-markovian chain that asymptotically converges to a flat-histogram Metropolis-Hastings. The Wang-Landau algorithm starts with an approximation of $G(E)$, $G_{WL}(E)=1/|E_{\max}-E_{\min}|$, and, on step 4 of the MH (see algorithm in Sec. III.2), it multiplies $G_{WL}(E(\boldsymbol{x}))$ ($\boldsymbol{x}$ is the current state of the random walk) by a constant $f>1$. After a given number of steps, $f$ is reduced by a factor 2 and this procedure is repeated until a final $f_{\min}\simeq 1$ is reached Wang and Landau (2001). The value $f_{\min}$ and how $f$ is reduced dictates how close $G_{WL}(E)$ will be from $G(E)$ Wang and Landau (2001); Zhou and Bhatt (2005); Belardinelli et al. (2008).

III.4 Characterisation of the efficiency

The relative error $\epsilon(E)$ depends on the value of $N$ and on the particular value of $E/N$. To avoid discussing the dependency on $E/N$, the efficiency of the flat-histogram is often quantified in terms of the average round-trip Trebst et al. (2004); Dayal et al. (2004); Viana Lopes (2006), which is an upper bound for the number of Markov steps $m$ (samples) required to obtain an uncorrelated sample from $\pi(\boldsymbol{x})$ Costa et al. (2005); Viana Lopes (2006). The round-trip, $\tau$, is the average number of steps required for the algorithm to go from a state $\boldsymbol{x}$ with $E(\boldsymbol{x})=E_{\min}$ to a state $\boldsymbol{x}$ with $E(\boldsymbol{x})=E_{\max}$ and return back. 333Numerically, the round-trip time is computed by having a boolean value tracking whether the random walk is moving in the direction $E_{\min}\rightarrow E_{\max}$ or $E_{\max}\rightarrow E_{\min}$, and use it to measure the total number of steps required to make the full round-trip. Ref. Leitão et al. (2013) shows how the autocorrelation time of a canonical ensemble is related to the autocorrelation time of a flat-histogram, and therefore we will use here the round-trip time of a flat-histogram simulation to quantify the efficiency for both distributions. The uniform sampling has a round-trip that increases exponentially with $N$: on average it takes $1/G(E_{\max})$ samples to get one sample with $E_{\max}$, and $1/G(E_{\max})$ increases exponentially with $N$.

Importance sampling Monte Carlo is widely used in statistical physics because the computational cost often scales polynomially with $N$ Viana Lopes (2006); Robert and Casella (2005), which is a dramatic improvement over uniform sampling. Under the hypothesis that (a) $\Delta E\equiv E(\boldsymbol{x}^{\prime})-E(\boldsymbol{x})\approx 1\ll N$ and (b) the correlation between the different $E(\boldsymbol{x})$ of the random walk decay fast enough, it can be shown that the roundtrip $\tau$ scales as Viana Lopes (2006); Leitão et al. (2014)

[TABLE]

For example, consider the problem of sampling states of an open chaotic system with different escape times $t_{e}$ from $t_{e}=1$ up to $t_{e}=t_{e}^{\max}$. In this case, $E_{\min}=1$ and $E_{\max}=N=t_{e}^{\max}$. Thus, under the hypothesis (a) and (b) above, the round-trip is expected to scale as

[TABLE]

There are known deviations of the scaling in Eq. 28 leading to a higher exponent, $N^{2+z}$ with $z>0$ Trebst et al. (2004); Dayal et al. (2004); Viana Lopes (2006); Leitão et al. (2014) and there are two common situations where this happens: (1) the acceptance rate decreases drastically with $N$ (hypothesis (a) above is violated); (2) autocorrelation drastically increases with $N$ (hypothesis (b) above is violated) Viana Lopes (2006); Leitão et al. (2014). Nevertheless, these do not qualitatively change the argument: Monte Carlo with flat-histogram has the potential of generating rare states polynomially with increasing $N$, while uniform sampling generates rare states exponentially with increasing $N$.

III.5 Summary: challenges for the application of Monte Carlo methods to chaotic systems

To achieve a MH algorithm that scales polynomially with N, the autocorrelation time needs to be low. A key ingredient for an efficient algorithm is a good proposal distribution $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ Roberts et al. (1997). The ideal proposal distribution of a Metropolis-Hastings algorithm draws $\boldsymbol{x}^{\prime}$ independently of $\boldsymbol{x}$ according to $\pi(\boldsymbol{x}^{\prime})$, $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})=\pi(\boldsymbol{x}^{\prime})$. This is because i) the acceptance in Eq. 23 is always one and ii) each step of the random walk generates an independent sample $\boldsymbol{x}$, which implies that the error in Eq. 22 is minimal. The difficulty of sampling rare events in chaotic systems is that a useful $\pi(\boldsymbol{x})$ to sample them is a difficult function to sample from. For concreteness, consider the problem of sampling high escape times in the open tent map defined in Appendix VII.1.3 – Eq. 66, whose $t_{e}(\boldsymbol{x})$ is represented in Fig. 12 – and consider the canonical sampling distribution, $\pi(\boldsymbol{x})\propto\exp(-\beta t_{e}(\boldsymbol{x}))$. In this example, $\pi(\boldsymbol{x})\propto\exp(-\beta 1)$ between $[1/a,1-1/a]$ and so forth ($\exp(-\beta t_{e})$) in subsequent intervals. The number of intervals increases as $2^{t_{e}}$. Therefore, sampling from $\pi(\boldsymbol{x})$ would require enumerating every interval, sample one at random according to a correct distribution that depends on $\beta$, and then sample a uniform point within that interval. While this can be done in simple maps such as the open tent map, this is unfeasible in a general system where the $t_{e}(\boldsymbol{x})$ dependency is unknown.

One way to approach the problem could be to consider $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ to be the uniform distribution over $\Gamma$. One could imagine that changing the sampling distribution $\pi(\boldsymbol{x})$ could decrease the variance of the estimator, since this gives more preference to rarer states. However, this is not the case: changing the sampling distribution alone does not decrease the scaling of the variance of the estimator. This is because changing the sampling distribution (e.g. using a canonical ensemble) leads to an exponential increase of the autocorrelation time $T(E)$ with $E$, see Appendix VII.2, making it as efficient as traditional uniform sampling.

In summary, Metropolis-Hastings is an excellent candidate to approach the numerical challenges found in the study of rare events in chaotic systems. Firstly, because it is grounded in strong mathematical results such as importance sampling theorem and asymptotic convergence of Markov processes. Secondly, because it is formulated with very little assumptions about the system, the observable of interest or the dynamics of the system, which gives enough freedom to adapt it to the specific aim (sampling or finding), observable, and system. Thirdly, because there seems to be no theoretical reason for the sampling to be exponential; Metropolis-Hastings is used to sample rare states in polynomial time in other problems of statistical physics. Finally, because the numerical problems found in chaotic systems can be re-written as problems where Metropolis-Hastings is suitable for. On the other hand, the optimal proposal of Metropolis-Hastings is unfeasible in chaotic systems, and without any extra information about the system, Metropolis-Hastings is as efficient as the traditional uniform sampling.

IV Proposal distribution

The problem we address in this section is: how to incorporate general properties of chaotic systems into the proposal distribution in such a way that the Metropolis-Hastings algorithm becomes efficient? We first set an aim for the proposal distribution and we then show how this aim can be achieved in the different problems involving deterministic chaotic system.

IV.1 Aim of the proposal distribution

The goal of the proposal distribution $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ we construct here will be to bound the acceptance rate in Eq. 23 away from 0 and 1. Since $a(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ in Eq. 23 depends on $\boldsymbol{x}^{\prime}$, it is essential to look at its expectation over the proposal,

[TABLE]

Our goal is to construct a proposal distribution such that444E.g. Ref. Roberts et al. (1997) computed $a_{\star}=0.234$.

[TABLE]

The motivation to set this as the starting point (and cornerstone) of our method is that it avoids the two typical origins of high correlation times $T$. When the acceptance is low, $T$ increases because the method remains stuck in the same state $\boldsymbol{x}$ for long times. High acceptance typically indicates that $\boldsymbol{x}^{\prime}$ is too close to $\boldsymbol{x}$ (often $E(\boldsymbol{x}^{\prime})=E(\boldsymbol{x})$), which implies that the simulation moves too slowly (in $E$ and in $\Omega)$. The goal here is not that the acceptance is exactly $a_{\star}$, but rather that it remains bounded from the extremes and that it does not strongly depend on $N$.

In general it is non-trivial to construct a proposal distribution that guarantees a constant acceptance. Thus, the next step is to approximate Eq. 31 by a simpler condition. Let us first notice that $\pi(\boldsymbol{x})$ only depends on $\boldsymbol{x}$ through $E_{\boldsymbol{x}}\equiv E(\boldsymbol{x})$, $\pi(\boldsymbol{x})=\pi(E_{\boldsymbol{x}})$. Therefore, at the very least, the proposal $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ should guarantee that $\boldsymbol{x}^{\prime}$ is generated from $\boldsymbol{x}$ in such a way that $\pi(E_{\boldsymbol{x}^{\prime}})$ is neither too close (high acceptance) nor too far (low acceptance) from $\pi(E_{\boldsymbol{x}})$. Quantitatively, this can be written as

[TABLE]

where $0<a<1$ is a constant. When the proposal is able to achieve a small variation of $E$, $\pi(E_{\boldsymbol{x}^{\prime}})$ can be expanded in Taylor series around $E_{\boldsymbol{x}^{\prime}}=E_{\boldsymbol{x}}$, which allows to write

[TABLE]

555In statistical physics an heuristics often used in Metropolis-Hastings is $\Delta E\equiv E_{\boldsymbol{x}^{\prime}}-E_{\boldsymbol{x}}\approx 1$. E.g. it is used in the derivation of Eq. 28, and has been used for example in spin systems (single spin flip) Viana Lopes (2006), ensemble of complex networks (single link exchange) Fischer2015, and proteins Grassberger1997.

Inserting Eq. 33 in Eq. 32, an average constant acceptance is thus achieved when

[TABLE]

This equation, the main result of this section, is a condition that an efficient Metropolis-Hastings imposes to the proposal distribution in terms of the average difference in the observable $E$. This condition is non-trivial because it depends on the particular $\pi$, $E$, $\boldsymbol{F}$, and $\boldsymbol{x}$. For the sampling distributions discussed in the previous section, we obtain:

Canonical ensemble:

when $\pi(\boldsymbol{x})=e^{-\beta E_{\boldsymbol{x}}}$, the condition in Eq. 34 is given by

[TABLE]

That is, the higher the $\beta$, the closer the proposed $E_{\boldsymbol{x}^{\prime}}$ has to be from $E_{\boldsymbol{x}}$.

Flat-histogram

When $\pi(\boldsymbol{x})\propto 1/G(E_{\boldsymbol{x}})$, the condition in Eq. 34 is given by

[TABLE]

When $E_{\boldsymbol{x}}$ is close to the maximum of $G(E)$, the derivative of $\log G$ approaches 0 and $E_{\boldsymbol{x}^{\prime}}$ can be arbitrary distant from $E_{\boldsymbol{x}}$. As $E_{\boldsymbol{x}}$ deviates from the maximum of $G(E)$, smaller changes are necessary to achieve a constant acceptance. Note that equation 34 is valid for the Metropolis-Hastings algorithm in general and should be of interest also in other contexts (e.g., arbitrarily large number of spins can be flipped close to the maximum of the density of states).

IV.2 Propose correlated trajectories

The proposal distribution requires correlating the trajectory starting at $\boldsymbol{x}^{\prime}$ with a trajectory starting at $\boldsymbol{x}$ such that Eq. 34 holds. Fulfilling this requirement requires the ability to control $\mathbb{E}\left[E_{\boldsymbol{x}^{\prime}}-E_{\boldsymbol{x}}|\boldsymbol{x}\right]$, which requires a procedure to correlate the state $\boldsymbol{x}^{\prime}$ with $\boldsymbol{x}$. The aim of this section is to introduce a quantification of the correlation of two trajectories of finite-time $t_{o}$ that can be related with $\mathbb{E}\left[E_{\boldsymbol{x}^{\prime}}-E_{\boldsymbol{x}}|\boldsymbol{x}\right]$, and present two different proposal distributions that propose a state $\boldsymbol{x}^{\prime}$ on which this correlation is controlled.

The observables $E$ introduced in section II.3 are all dependent not only of $\boldsymbol{x}$, but also of the full trajectory of length $t_{o}$ starting $\boldsymbol{x}$.666The escape time $t_{e}$ requires that $\boldsymbol{F}^{t}(\boldsymbol{x})\notin\Lambda\ \ ,\forall t$; the FTLE $\lambda_{t_{o}}$ is a sum of terms computed over the trajectory. One natural way to quantify the similarity of two trajectories is the length the two trajectories remain within a distance $\Delta$ much smaller than a characteristic length of $\Gamma$. Formally, this can be quantified by $t_{\star}(\boldsymbol{x},\boldsymbol{x}^{\prime})=\max\{t\leq t_{o}:|\boldsymbol{F}^{t}(\boldsymbol{x})-\boldsymbol{F}^{t}(\boldsymbol{x}^{\prime})|\leq\Delta\}$. Under this definition, $0\leq t_{\star}(\boldsymbol{x},\boldsymbol{x}^{\prime})\leq t_{o}$. When $\boldsymbol{x}^{\prime}=\boldsymbol{x}$, $t_{\star}(\boldsymbol{x},\boldsymbol{x}^{\prime})=t_{o}$ because the trajectories are the same; when $\boldsymbol{x}^{\prime}$ is far from $\boldsymbol{x}$, $t_{\star}(\boldsymbol{x},\boldsymbol{x}^{\prime})=0$. However, there is another possibility for two trajectories starting at $\boldsymbol{x}^{\prime}$ and $\boldsymbol{x}$ to be similar: when the two trajectories are similar apart from a shift in time, see right panel of figure 3. One way to include both cases in our measure of similarity is to define $t_{\star}(\boldsymbol{x},\boldsymbol{x}^{\prime})$ as

[TABLE]

For ${t_{\text{shift}}}=0$, this recovers the case of two trajectories starting close to each other; ${t_{\text{shift}}}\neq 0$, it includes situations where a trajectory starts at $\boldsymbol{x}^{\prime}$ close to $\boldsymbol{F}^{{t_{\text{shift}}}}(\boldsymbol{x})$.

This definition is motivated by the concept of symbolic sequences. Ott (1993) The similarity of the trajectory starting from $\boldsymbol{x}^{\prime}$ with the one starting from $\boldsymbol{x}$ can be quantified by the number of symbols that both trajectories share, which corresponds to the $t_{\star}(\boldsymbol{x},\boldsymbol{x}^{\prime})$ in Eq. 37. The definition in Eq. 37 avoids the necessity of the existence of a phase-space partition, but, for the purposes of the argument below, the two trajectories share a sequence of $t_{\star}$ states that are close within $\Delta$.

The average correlation between two states whose one is drawn according to a proposal distribution is here defined by

[TABLE]

Notice that this quantity does not depend on the particular sampling distribution or algorithm; it is a function of the proposal distribution and the state $\boldsymbol{x}$. The goal of the next sub-sections is to construct proposal distributions that guarantee a given average correlation $t_{\star}(\boldsymbol{x})$. We will show, for example, that we can enforce an average correlation $t_{\star}(\boldsymbol{x})$ if we use a normal distribution centered around $\boldsymbol{x}$ with a specific standard deviation as our proposal distribution.

IV.2.1 Shift proposal

One proposal that guarantees that trajectories are correlated by $t_{\star}$ is the shift proposal, originally introduced in Ref. Dellago et al. (2002) in the context of sampling paths of chemical reactions, and has also been used in Ref Grünwald et al. (2008). It consists in proposing a state $\boldsymbol{x}^{\prime}$ that is a forward or backward iteration of $\boldsymbol{x}$, $\boldsymbol{x}^{\prime}=\boldsymbol{F}^{t_{\text{shift}}}(\boldsymbol{x})$, where ${t_{\text{shift}}}$ is a free parameter. The relation between ${t_{\text{shift}}}={t_{\text{shift}}}(\boldsymbol{x})$ and $t_{\star}=t_{\star}(\boldsymbol{x})$ is that a shift of $\pm{t_{\text{shift}}}$ guarantees that $t_{o}-{t_{\text{shift}}}$ elements of the original trajectory are preserved. Therefore, this proposal guarantees that $t_{\star}$ elements are preserved when $|{t_{\text{shift}}}|=t_{o}-t_{\star}$, see right panel of Fig. 3. Because detailed balance has to be enforceable, the proposal must contain backward and forward shifts. A proposal that automatically fulfils detailed balance is one on which the backward and forward shifts are equally likely:

[TABLE]

with ${t_{\text{shift}}}={t_{\text{shift}}}(\boldsymbol{x})=t_{o}-t_{\star}(\boldsymbol{x})$. Given the target average correlation $t_{\star}(\boldsymbol{x})$, this proposal can be implemented as follows: generate a random number $r\in[0,1]$; if $r<0.5$, make $\boldsymbol{x}^{\prime}=F^{t_{o}-t_{\star}(\boldsymbol{x})}(\boldsymbol{x})$, else, make $\boldsymbol{x}^{\prime}=F^{t_{\star}(\boldsymbol{x})-t_{o}}(\boldsymbol{x})$.

This proposal unfortunately has some disadvantages: i) a priori there is no guarantee that $\boldsymbol{F}^{t_{\text{shift}}}(\boldsymbol{x})\in\Gamma$. It is applicable when $\Gamma=\Omega$, which e.g. is not the case in open systems; ii) it requires the map to be invertible; iii) the proposal can only propose states that are forward or backward iterations of $\boldsymbol{x}$. Consequently, for the random walk to be ergodic in the phase-space, the map itself must be ergodic. iv) this proposal diffuses without drift on a trajectory passing through $\boldsymbol{x}$, by shifting the starting point forward or backward. Thus, it will always sample fewer states than a time average of a trajectory starting at $\boldsymbol{x}$. On the other hand, the main advantage of this proposal is that it performs non-local jumps in the phase-space. That is, it allows to jump from a region of the phase-space to another region while still maintaining $\boldsymbol{x}^{\prime}$ correlated with $\boldsymbol{x}$. As shown below, in combination with other proposal, this proposal is useful to reduce correlations stemmed from local jumps.

IV.2.2 Neighbourhood proposal

Another strategy to construct a proposal on which on average the states are correlated by $t_{\star}(\boldsymbol{x})$ is to perturb $\boldsymbol{x}$ by a finite amount $\boldsymbol{\delta}$, $\boldsymbol{x}^{\prime}=\boldsymbol{x}+\boldsymbol{\delta}$, characterised by a direction $\hat{\boldsymbol{\delta}}$ and a norm $\delta$, $\boldsymbol{\delta}\equiv\hat{\boldsymbol{\delta}}\delta$. A common case is when the probability distribution is separated in two independent terms Robert and Casella (2005):

[TABLE]

and that $P(\hat{\boldsymbol{\delta}}|\boldsymbol{x})$ is uniformly distributed in the $D$ directions and $P(\delta|\boldsymbol{x})$ has zero mean (i.e. an isotropic proposal). Here we restrict the analysis to this situation, and we also assume that $P(\delta|\boldsymbol{x})$ is characterised by a well defined scale, e.g. it is an half-normal distribution777An exponential distribution would be equally acceptable and would not change the main conclusions. with mean $\delta_{x}(\boldsymbol{x})$:

[TABLE]

This choice makes the ratio $g(\boldsymbol{x}|\boldsymbol{x}^{\prime})/g(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ in Eq. 23 to be

[TABLE]

The main motivation for this choice is that the proposal distribution is described by a single function, $\delta_{x}(\boldsymbol{x})$, that quantifies the distance $\boldsymbol{x}^{\prime}-\boldsymbol{x}$, $\mathbb{E}\left[|\boldsymbol{x}^{\prime}-\boldsymbol{x}||\boldsymbol{x}\right]=\delta_{x}(\boldsymbol{x})$.

The goal is now to relate $\delta_{x}(\boldsymbol{x})$ with $t_{\star}(\boldsymbol{x})$. Let us start to describe two important limits: in the limit $\delta_{x}(\boldsymbol{x})\rightarrow 0$, the states are the same and therefore $\lim_{\delta_{x}(\boldsymbol{x})\rightarrow 0}t_{\star}(\boldsymbol{x})=t_{o}$. In the limit $\delta_{x}(\boldsymbol{x})\rightarrow|\Gamma|$, the proposal is approximately equal to draw $\boldsymbol{x}^{\prime}$ uniformly from $\Gamma$, and $\boldsymbol{x}^{\prime}$ is independent of $\boldsymbol{x}$ and $t_{\star}(\boldsymbol{x})=0$. To preserve a correlation of $t_{\star}(\boldsymbol{x})$, it is necessary that $\delta_{x}(\boldsymbol{x})$ is such that the two trajectories starting at $\boldsymbol{x}$ and $\boldsymbol{x}^{\prime}$ are close together up a time $t_{\star}(\boldsymbol{x})$, see left panel of Fig. 3. Because the system is chaotic, for small $\delta_{x}(\boldsymbol{x})$, two trajectories diverge exponentially in time according to Eq. 4, and, in particular, their maximal distance is given by Eq. 7. Therefore, to guarantee that two trajectories are distanced at most by $\Delta$ after a time $t_{\star}(\boldsymbol{x})$, $\delta_{x}(\boldsymbol{x})$ must be given by

[TABLE]

This equation relates the parameter of the proposal distribution, $\delta_{x}(\boldsymbol{x})$, with the average correlation $t_{\star}(\boldsymbol{x})$ of the two states $\boldsymbol{x}$ and $\boldsymbol{x}^{\prime}$.

The neighbourhood proposal derived above is closely related to a proposal described in Ref. Grünwald et al. (2008) as “precision shooting”. Precision shooting constructs a trajectory $\{\boldsymbol{x}_{i}^{\prime}\}$ with $\boldsymbol{x}_{0}^{\prime}\equiv\boldsymbol{x}^{\prime}=\boldsymbol{x}+\delta\hat{\boldsymbol{\delta}}$ (where $\delta$ is a free parameter) that, within the numerical precision of a computer, is indistinguishable (in the system considered) from a trajectory starting at $\boldsymbol{x}^{\prime}$ with $\delta$ small. The trajectory $\{\boldsymbol{x}_{i}^{\prime}\}$ shadows a true trajectory starting at $\boldsymbol{x}^{\prime}$, in the same spirit as the algorithm used in Ref. Grebogi et al. (1990) to construct a pseudo-trajectory. Thus, precision shooting can be interpreted as the neighbourhood proposal, Eq. 43, with $t_{\star}$ a free parameter (related to $\delta$ via Eq. 43), that assumes shadowing theorem to simplify the construction $\boldsymbol{x}^{\prime}$. Ref. Grünwald et al. (2008) discusses how the acceptance rate depends on $\delta$, suggesting that the acceptance rate increases with decreasing $\delta$ (Fig. 9 of the ref.). In light of the discussion in section IV.2, this result is interpreted as follows: as $\delta$ decreases, $\boldsymbol{x}^{\prime}$ becomes more correlated with $\boldsymbol{x}$ (since $t_{\star}$ is related with $\delta$ by Eq. 43), and therefore the acceptance is expected to increase, as indicated in Fig. 9 of Ref. Grünwald et al. (2008). This discussion is unfortunately insufficient to us because it does not allow to derive $\delta$ (or $t_{\star}$) that fulfils the condition in Eq. 34. The crucial advantage of Eq. 43 is that it allows to relate $\mathbb{E}\left[E(\boldsymbol{x}^{\prime})-E(\boldsymbol{x})|\boldsymbol{x}\right]$ (in Eq. 34) with $\boldsymbol{x}^{\prime}-\boldsymbol{x}$ (in Eq. 43) through $t_{\star}$ (in Eq. 37). This is the goal of the next section.

IV.3 Guarantee local proposals

Now that we derived proposal distributions that enforce an average correlation $t_{\star}(\boldsymbol{x})$, the next (and final) step is to obtain a relationship between $t_{\star}(\boldsymbol{x})$ and $\mathbb{E}\left[E(\boldsymbol{x}^{\prime})-E(\boldsymbol{x})|\boldsymbol{x}\right]$. Because the computation of $t_{\star}(\boldsymbol{x})$ depends on the particular observable $E$, a different derivation is presented for two observables, $t_{e}$ and $\lambda_{t_{o}}$. Given the limitations of the shift proposal described above, the argumentation below is limited to neighbourhood proposals (an equivalent argumentation can be made to the shift proposal).

IV.3.1 FTLE in closed systems

As introduced in section II.1, the observable in this case is given by $E(\boldsymbol{x})=t_{o}\lambda_{t_{o}}(\boldsymbol{x})$, where $t_{o}$ is the finite time. The aim in this case is thus to write $\mathbb{E}\left[t_{o}\lambda_{t_{o}}(\boldsymbol{x}^{\prime})-t_{o}\lambda_{t_{o}}(\boldsymbol{x})|\boldsymbol{x}\right]$ as a function of $t_{\star}(\boldsymbol{x})$. The finite-time Lyapunov exponent considered in Eq. 8 is a sum of $t_{o}$ terms888This equation is strictly valid for 1D, but the idea is that we split the FTLE between a component up to $t_{\star}$ and the remaining component. and thus $t_{o}\lambda_{t_{o}}(\boldsymbol{x})$ can be written as the sum of the FTLE up to time $t_{\star}$ and the FTLE from $t_{\star}$ up to $t_{o}$,

[TABLE]

where $\boldsymbol{x}_{t_{\star}}\equiv\boldsymbol{F}^{t_{\star}}(\boldsymbol{x})$. Likewise for the the trajectory starting from $\boldsymbol{x}^{\prime}$,

[TABLE]

Because $\boldsymbol{x}^{\prime}$ is proposed according to Eq. 43, by construction, the first $t_{\star}$ states of the trajectory starting at $\boldsymbol{x}^{\prime}$ are close (within $\Delta$) to the states of the trajectory starting at $\boldsymbol{x}$ up to $t_{\star}(\boldsymbol{x})$. Therefore, we can approximate that the respective Lyapunovs up to time $t_{\star}$ are equal,

[TABLE]

Subtracting Eq. 44 from Eq. 45 and using Eq. 46 gives

[TABLE]

The left side of this equation is the same as in Eq. 34 and thus the aim now is to write the right side as a function of properties of the system. Let us focus on the calculation of $\mathbb{E}\left[\lambda_{t_{o}-t_{\star}}(\boldsymbol{x}_{t_{\star}}^{\prime})|\boldsymbol{x}\right]$ first. By construction, $\boldsymbol{x}^{\prime}$ is generated such that $|\boldsymbol{x}_{t_{\star}}^{\prime}-\boldsymbol{x}_{t_{\star}}|\approx\Delta$. Because the system is chaotic, one can approximate that $\boldsymbol{x}_{t_{\star}}^{\prime}$ is sufficiently separated from $\boldsymbol{x}_{t_{\star}}$ such that $\lambda_{t_{o}-t_{\star}}(\boldsymbol{x}_{t_{\star}}^{\prime})$ is independent of $\boldsymbol{x}$. Under this approximation, $\boldsymbol{x}_{t_{\star}}^{\prime}$ is essentially a random state from the phase-space, and thus $\lambda_{t_{o}-t_{\star}}(\boldsymbol{x}_{t_{\star}}^{\prime})$ will be a drawn from the distribution of FTLE at time $t_{o}-t_{\star}$. Denoting the mean of this distribution by $\lambda_{L,t_{o}-t_{\star}}$, we get

[TABLE]

Replacing Eq. 48 in Eq. 47 gives

[TABLE]

This equation relates the expected change in the observable with properties of the system ($\lambda_{L,t_{o}-t_{\star}}$), of the trajectory $\boldsymbol{x}$, $\lambda_{t_{o}-t_{\star}}(\boldsymbol{x}_{t_{\star}})$, and $t_{\star}(\boldsymbol{x})$ and it can thus be used in the energy condition we obtained earlier, Eq. 34. Replacing the left side of Eq. 34 by the expectation in Eq. 49 and solving to $t_{\star}(\boldsymbol{x})$ gives

[TABLE]

We can further simplify this relation with two approximations: a) for large $t_{o}-t_{\star}$, the mean FTLE at time $t_{o}-t_{\star}$, $\lambda_{L,t_{o}-t_{\star}}$, is approximately the Lyapunov exponent of the system, $\lambda_{L}$,

[TABLE]

b) because trajectories of chaotic systems are short-correlated in time, the FTLE of the trajectory $\boldsymbol{x}$ up to $t_{\star}$ will be approximately equal to the FTLE of the trajectory starting at $\boldsymbol{x}_{t_{\star}}$. Thus,

[TABLE]

Using Eq. 51 and Eq. 52, $t_{\star}(\boldsymbol{x})$ in Eq. 50 can be written as

[TABLE]

To enforce that $t_{\star}\in[0,t_{o}]$, we use

[TABLE]

which is the main result of this section. This equation provides an expression to $t_{\star}(\boldsymbol{x})$ that can be inserted in the parameter of the proposal distribution, Eq. 43, that under the approximations used above achieves a constant acceptance rate999 In Ref. Leitão et al. (2014) we used $t_{\star}=t_{o}-1$, different from Eq. 54. The derivation of Eq. 54 uses stronger approximations than the derivation of $t_{\star}=t_{o}-1$. See argumentation after Eq. (14) of Ref. Leitão et al. (2014). .

The derivation of Eq. 54 can be generalised to other observables which, as $\lambda_{t_{o}}$, can be written as an average over the trajectory: consider

[TABLE]

where $f(\boldsymbol{x}_{i})$ is an arbitrary function of the phase-space (the logarithm of the derivative of the map corresponds to $e_{t_{o}}(\boldsymbol{x})=\lambda_{t_{o}}(\boldsymbol{x})$, $E_{t_{o}}(\boldsymbol{x})=\lambda_{t_{o}}(\boldsymbol{x})t_{o}$). Replacing this quantity in the derivation of Eq. 54 mutatis mutandis and without using the approximation in Eq. 52, one obtains

[TABLE]

where $E_{L}$ is approximately the mean of the distribution $P(E_{t_{o}})$ (using the approximation in Eq. 48). This generalizes Eq. 54 for an arbitrary average over trajectories of size $t_{o}$, $e_{t_{o}}(\boldsymbol{x})$, and it should be useful to sample rare states in respect to observables correspondent to expected values over trajectories.

IV.3.2 Escape time in strongly chaotic open systems

As introduced before, in strongly chaotic open systems $E(\boldsymbol{x})=t_{e}(\boldsymbol{x})$ and $P(t_{e})\sim\exp(-\kappa t_{e})$. The aim here is to compute $t_{\star}(\boldsymbol{x})$ that fulfils Eq. 34 by taking into account that $\boldsymbol{x}$ is given and $\boldsymbol{x}^{\prime}=\boldsymbol{x}+\hat{\boldsymbol{\delta}}\delta_{x}(\boldsymbol{x})$ with $\delta_{x}(\boldsymbol{x})$ given by Eq. 43. In this case the trajectory’s length is not a constant but instead it is given by the escape time $t_{e}(\boldsymbol{x})$.

A trajectory starting at $\boldsymbol{x}^{\prime}$ proposed according to Eq. 43 fulfils $|\boldsymbol{x}_{t_{\star}}^{\prime}-\boldsymbol{x}_{t_{\star}}|\leq\Delta$. Therefore, up to $t_{\star}$, the two trajectories are indistinguishable. Under this assumption, and because from the definition of $t_{\star}$ in Eq. 37, $t_{e}(\boldsymbol{x}^{\prime})=t_{\star}(\boldsymbol{x})+t_{e}(\boldsymbol{x}_{t_{\star}}^{\prime})$ and therefore

[TABLE]

where $t_{e}(\boldsymbol{x}_{t_{\star}}^{\prime})$ is the escape time of the state $\boldsymbol{x}^{\prime}$ iterated $t_{\star}$ times and $\mathbb{E}\left[t_{e}(\boldsymbol{x}_{t_{\star}}^{\prime})|\boldsymbol{x}\right]=\int_{\Omega}g(\boldsymbol{x}^{\prime}|\boldsymbol{x})t_{e}(\boldsymbol{x}_{t_{\star}}^{\prime}){\bf d}\boldsymbol{x}^{\prime}$ is the expected $t_{e}(\boldsymbol{x}_{t_{\star}}^{\prime})$ over a neighbourhood of $\boldsymbol{x}$ of size $\delta_{x}(\boldsymbol{x})$ (given by Eq. 43)).

The idea behind this equation is represented in Fig. 4. The proposal density $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ at $\boldsymbol{x}$ is such that, at $t_{\star}$, the two trajectories distance themselves on average by $\Delta$. After $t_{\star}$, the trajectory starting at $\boldsymbol{x}^{\prime}_{t_{\star}}$ will approximately be independent of $\boldsymbol{x}$ and thus $t_{e}(\boldsymbol{x}^{\prime})=t_{\star}+t_{e}(\boldsymbol{x}^{\prime}_{t_{\star}})$. Moreover, proposing according to Eq. 43 is equivalent to “zoom” the landscape of $t_{e}(\boldsymbol{x})$ around $\boldsymbol{x}$ with a scale correspondent to $t_{\star}$’th iteration of the construction of the landscape. Given the self-similarity of the landscape (see e.g. Fig. 13), under this zoom, the landscape of $t_{e}(\boldsymbol{x}_{t_{\star}})$ is equal to the landscape of $t_{e}(\boldsymbol{x})$ and therefore $\mathbb{E}\left[t_{e}(\boldsymbol{x}_{t_{\star}}^{\prime})|\boldsymbol{x}\right]$ should be a constant independent of $t_{e}$. Moreover, because $\boldsymbol{x}_{t_{\star}}^{\prime}$ is approximately independent of $\boldsymbol{x}$, $\mathbb{E}\left[t_{e}(\boldsymbol{x}_{t_{\star}}^{\prime})|\boldsymbol{x}\right]$ is the just average escape time of an independent state $\boldsymbol{x}^{\prime}_{t_{\star}}$, which is the average of $P(t_{e})$ and is given approximately by $1/\kappa$. Thus, $\mathbb{E}\left[t_{e}(\boldsymbol{x}_{t_{\star}}^{\prime})|\boldsymbol{x}\right]$ is the average escape time of an independent state $\boldsymbol{x}^{\prime}_{t_{\star}}$, which is the average of $P(t_{e})$ and is given by $1/\kappa$:

[TABLE]

It is this result that incorporates the self-similarity of the escape time function: the value of $t_{\star}$ chooses the particular zoom of the landscape (Fig. 4), and this equation assumes that, as long as the zoom is proportional to $\exp(\lambda t_{\star})$, the average escape time of the phase-space of the zoomed region is still $1/\kappa$. Replacing Eq. 58 in Eq. 57 gives

[TABLE]

and subtracting $t_{e}(\boldsymbol{x})$ on both sides of Eq. 59 gives

[TABLE]

The left side of this equation is the left side of the condition of constant acceptance rate, Eq. 34. Equating both left sides and solving for $t_{\star}(\boldsymbol{x})$ gives

[TABLE]

which is analogous to Eq. 54 and is the central result of this section. Together with Eq. 43, it is the proposal distribution we had the goal of constructing.

IV.4 Summary: how to propose

The argument in this section can be summarised as follows:

• •

Metropolis-Hastings requires a proposal that guarantees a specific variation of $E$, which we approximate by Eq. 34.

• •

the proposal with a scale given by Eq. 43 guarantees that on average the trajectory starting at $\boldsymbol{x}^{\prime}$ stays close to the trajectory starting at $\boldsymbol{x}$ up to $t_{\star}(\boldsymbol{x})$.

• •

the variation of $E$ is related to $t_{\star}(\boldsymbol{x})$ via Eqs. 49 and 60.

These three steps led to analytical formulas for $t_{\star}(\boldsymbol{x})$, Eqs. 61 and 54, that make the proposal satisfy Eq. 34. Together with Eq. 43, they define proposal distributions required for an average constant acceptance.

There are 3 points that deserve to be noted: the first point is that the formulas we obtained for $t_{\star}(\boldsymbol{x})$ are also valid for the problem of finding rare states discussed e.g. in Refs. Sweet et al. (2001); Tailleur and Kurchan (2007); Kitajima and Iba (2011); Iba et al. (2014); Bollt (2005). Specifically, these formulas were derived to guarantee Eq. 34, which dictates how different $\mathbb{E}\left[E(\boldsymbol{x}^{\prime})|\boldsymbol{x}\right]$ has to be from $E(\boldsymbol{x})$ to guarantee a constant acceptance. This condition is stronger than the condition required for an algorithm to find minima or maxima of $E$, which requires only proposing states $\boldsymbol{x}^{\prime}$ such that $\mathbb{E}\left[E(\boldsymbol{x}^{\prime})|\boldsymbol{x}\right]>E(\boldsymbol{x})$ (or vice-versa for minimising $E$). This is independent of the particular minimisation algorithm (e.g. stimulated annealing, step descent, stagger and dagger in open systems) because it only discusses which new state $\boldsymbol{x}^{\prime}$ should be tried, given the current state $\boldsymbol{x}$. The second point is that Eqs. 61 and 54 reduce the proposal to a uniform distribution when the sampling distribution is the uniform distribution: $d\log\pi(E)/dE(E)\rightarrow 0$ implies $\delta_{x}(\boldsymbol{x})\rightarrow\infty$. The third point is that different proposals in the literature, precision shooting Grünwald et al. (2008); Leitão et al. (2014), exponential proposal distribution Sweet et al. (2001); Kitajima and Iba (2011), and the one in Ref. Leitão et al. (2013), can be obtained from the proposals derived in this section through different approximations, as shown in Appendix VII.3.

Overall, this section described a framework to add information of chaotic systems, in this case the self-similar properties of the landscape, the exponential divergence of trajectories, and the exponential decay of correlations, to construct an efficient Metropolis-Hastings algorithm to sample them. A simplified description of the algoirthm is given in in Appendix VII.4 and an open-source implementation in Ref. Leitão (2017). The next section is devoted to test the assumptions used here on each of the problems, escape time and FTLE, and confirm the practical usefulness of the framework.

V Numerical Tests

The previous section concluded with a set of formulas – the proposal in Eq. 43 combined with the formulas for $t_{\star}(\boldsymbol{x})$ Eq. 54 or Eq. 61 – for proposing states $\boldsymbol{x}^{\prime}$ that are expected to guarantee the desired acceptance rate (bounded from [math] and $1$). In this section we test some of the approximations made in the derivation of $t_{\star}(\boldsymbol{x})$ (in a simple system) and we analyse the efficiency of the algorithm (confirming the polynomial scaling) in the computation of the FTLE in closed system – introduced in Sec. II.1 – and of the escape time in open systems – introduced in Sec. II.2. The tests are performed in the skewed tent map and open skewed tent map, and the efficiency is tested in numerous maps (chain of couple Hénon maps, Standard map, Logistic map). All these systems are introduced in detail in Appendix VII.1.

V.1 Finite-time Lyapunov exponent

The first approximation made in the derivation of Eq. 54 is that when $\delta_{x}\equiv|\boldsymbol{x}^{\prime}-\boldsymbol{x}|$ is drawn from a half-normal distribution with scale parameter $\delta_{x}(\boldsymbol{x})$ given by Eq. 43, $\boldsymbol{F}^{t_{\star}}(\boldsymbol{x}^{\prime})$ is sufficiently close from $\boldsymbol{F}^{t_{\star}}(\boldsymbol{x})$ such that $\mathbb{E}\left[\lambda_{t_{\star}}(\boldsymbol{x}^{\prime})|\boldsymbol{x}\right]=\lambda_{t_{\star}}(\boldsymbol{x})$, Eq. 46, holds. We test this numerically by randomly drawing $10^{5}$ states $\boldsymbol{x}_{i}$ in the tent map with $a=3$ (see Appendix VII.1), and, for each, propose a state $\boldsymbol{x}_{i}^{\prime}=\boldsymbol{x}_{i}+\boldsymbol{\delta}(\boldsymbol{x})$ according to Eq. 41 with $\delta_{x}(\boldsymbol{x})$ given by our Eq. 43. From the $10^{5}$ pairs of states $(\boldsymbol{x}_{i},\boldsymbol{x}^{\prime}_{i})$, we estimate the difference of the observables: $E^{\prime}\equiv\mathbb{E}\left[t\lambda_{t}(\boldsymbol{x}_{i}^{\prime})|\boldsymbol{x}\right]$ and $E\equiv t\lambda_{t}(\boldsymbol{x}_{i})$. Our expectation is that for a fixed $t_{\star}=t$, $E^{\prime}-E$ should be much smaller than $E$ (the relevant scale in Eq. 47). We numerically obtain that within the $99\%$ quantile, $E^{\prime}-E\approx 1$ independently of $E$ and $t$. This value is much smaller than $E\in[0.6t,1.1t]$, specially since we are interested in large $t$. This strongly supports the approximation we make in Eq. 46.

The second approximation made in the derivation of Eq. 54 is that there is no dependence between $\lambda_{t_{\star}}(\boldsymbol{x})$ and $\lambda_{t_{o}-t_{\star}}(\boldsymbol{F}^{t_{\star}}(\boldsymbol{x}))$, Eq. 48. In other words, that the FTLE of the trajectory starting at $\boldsymbol{F}^{t_{\star}}(\boldsymbol{x})$ and ending at $t_{o}$ is indistinguishable from the one drawn from the distribution of FTLE with finite-time $t_{o}-t_{\star}$, $P(\lambda_{t_{o}-t_{\star}})$. This approximation was numerically tested by drawing points $\boldsymbol{x}_{i}$, computing the pairs $\lambda_{t}(\boldsymbol{x}_{i}),\lambda_{t}(\boldsymbol{F}^{t}(\boldsymbol{x}_{i}))$ (i.e. $2t_{\star}=t_{o}=2t$), and testing whether the conditional probability of $\lambda_{t}(\boldsymbol{F}^{t}(\boldsymbol{x}_{i}))$ equals the (unconditional) probability of $\lambda_{t}(\boldsymbol{x}_{i})$. The results in Fig. 5 confirm the equality of these probabilities.

We finally test whether a proposal using $t_{\star}(\boldsymbol{x})$ given by Eq. 54 guarantees a constant acceptance rate, the original motivation for our calculation. The test consisted in sampling $10^{6}$ states according to the following procedure: 1) uniformly draw a state $\boldsymbol{x}_{i}$ and compute $\lambda_{i}\equiv\lambda_{t_{o}}(\boldsymbol{x}_{i})$; 2) generate a state $\boldsymbol{x}^{\prime}_{i}$ according to the proposal distribution Eq. 41 with $\delta_{x}(\boldsymbol{x})$ given by Eq. 43 and $t_{\star}$ given by Eq. 54 ($\delta_{0}=1$), and compute $\lambda_{i}^{\prime}\equiv\lambda_{t_{o}}(\boldsymbol{x}^{\prime}_{i})$; 3) store $G_{i}\equiv g(\boldsymbol{x}|\boldsymbol{x}^{\prime})/g(\boldsymbol{x}^{\prime}|\boldsymbol{x})$ computed from Eq. 42: $\delta_{x}$ is given by Eq. 43, and $|\boldsymbol{x}^{\prime}-\boldsymbol{x}|$ is given by storing $\delta_{i}=|\boldsymbol{x}^{\prime}_{i}-\boldsymbol{x}_{i}|$. The ratio of the target distribution is given by $r_{i}\equiv\pi(E^{\prime})/\pi(E)=\exp(-\beta t_{o}(\lambda_{i}^{\prime}-\lambda_{i}))$ for the canonical ensemble and $r_{i}=G(E)/G(E^{\prime})$, where $G(E)=G(t_{o}\lambda_{i})$, is given by Eq. 64 for the flat-histogram. The numerically estimated acceptance ratio, $A(\lambda_{t})\equiv\left\langle\min(1,r_{i}\times G_{i})\right\rangle$ is shown in Fig. 6. It is not independent of $\lambda_{t_{o}}$ – the expected outcome based on the assumption of constant acceptance used in our derivations – but it is bounded from $0,1$ for increasing $N=t_{o}$ – the original requirement for an efficient simulation set in Sec. IV.1. In the canonical ensemble, there is linear dependency of $\Pi$ with $\lambda_{t_{o}}$, and in the flat-histogram ensemble, the ratio is 0.8 in the maximum of $P(\lambda_{t_{o}})$, and decays to about $0.1$ on the tails.

V.1.1 Efficiency of the flat-histogram ensemble

The success in achieving a bounded acceptance independent of $N$, as the one confirmed in the previous section for the tent map, does not guarantee the efficiency of the method (see Sec. III.4). Here we test the efficiency of flat-histogram computations of $P(\lambda_{t})$ in the tent map that use the neighbourhood proposal, the shift proposal, and both (mixed proposal). The results shown in Fig. 7 suggest that only when both proposals are used the efficiency scales polynomially with $N=t_{o}$.

This result can be understood looking at the landscape of $\lambda_{t}$, as illustrated in Figure 8. Imagine a flat-histogram simulation on this system, for $t_{o}=4$ (black curve in the figure), and analyse what happens to it in terms of a round-trip. Lets suppose that the simulation was recently at the minimum $\lambda_{t}$ (0.41) and that the next round-trip is made by going to the maximum $\lambda_{t}$ (1.09) and return back. Lets further suppose that the simulation eventually got to a state with $\lambda_{t}\approx 0.92$. Because $\pi(\boldsymbol{x})=\pi(\lambda_{t}(\boldsymbol{x}))$, every state at that $\lambda_{t}$ is equiprobable. Therefore, the state can be at any plateau (of the 4, see fig.), proportionally to their plateau-size. However, not every plateau contains, on its neighbourhood, a neighbour plateau with higher $\lambda_{t}$, for example, the plateau around $0.3$. Therefore, a local proposal would never be able to reach a higher plateau from a state on such a plateau. First, it would need to go backward, reach the maximum of $P(\lambda_{t})$ (around $0.69$), where the proposal proposes any other state, and then try to find another path towards a higher $\lambda_{t}$. This would already be the problem if the simulation would be a canonical ensemble with a $\beta$ favouring higher $\lambda_{t}$’s, since it would require decreasing $\lambda_{t}$ by an amount $\Delta\lambda$, and this happens with a probability that decays exponentially with $\Delta\lambda$. This is solved by using a flat-histogram ensemble. However, the crucial challenge here is that as $t$ increases, the number of local maxima also increases, but the number of maxima connected with the global maximum is constant: the red curve, with $t=6$, contains now $11$ plateaus for $\lambda_{t}\approx 0.87$, but only $1$ is locally connected to the maximum $\lambda_{t}$, the one around 0. This means that, as $t$ increases, it becomes more difficult to perform a round-trip: not only because the expected time to diffuse increases (see Eq. 28), but also because there are more times where the simulation diffuses forward and backward until it reaches the global maxima. The shift proposal alleviates this problem by allowing non-local proposals in the phase-space. A shift proposes a state $\boldsymbol{x}^{\prime}$ on a non-neighbourhood of $\boldsymbol{x}$, which improves the probability of reaching higher or lower $\lambda$’s, which explains why a combined proposal has such a low round-trip time in this system.

Finally, we confirm that a polynomial efficiency is obtained in the computation of $P(\lambda_{t_{o}})$ more generally. We performed flat-histogram simulations, using the Wang-Landau algorithm to estimate $P(\lambda_{t_{o}})$ on different chaotic systems: the tent map, the logistic map, and the standard map (see appendix VII.1 for details). The results are shown in Fig. 9 and confirm the dramatic improvement and generality of using Metropolis-Hastings to sample rare states. The simulations in Fig. 9 use $t_{\star}(\boldsymbol{x})=t_{o}-1$ instead of the one given by Eq. 54. This is computationally always more expensive because the correlations due to the neighbourhood proposal are maximal (see the discussion after Eq. 54), but on the other hand this proposal is simpler because it does not require estimating $d\log P/dE$ and $\lambda_{L}$.

V.2 Transient chaos

The derivation of Eq. 61 uses the assumption that when $\boldsymbol{x}^{\prime}$ is proposed with a scale given by Eq. 43 with $t_{\star}=t_{e}$, $\mathbb{E}\left[t_{e}(\boldsymbol{x}^{\prime})|\boldsymbol{x}\right]=t_{e}(\boldsymbol{x})$, as per Eq. 57. We tested assumption in a similar way we tested the approximation of Eq. 46 for the FTLE, and consisted in uniformly drawing states $\boldsymbol{x}_{i}$, compute their escape time $t_{e}\equiv t_{e}(\boldsymbol{x}_{i})$ and, for each, generate a state $\boldsymbol{x}_{i}^{\prime}=\boldsymbol{x}_{i}+\boldsymbol{\delta}(\boldsymbol{x})$ with a scale $\delta_{x}(\boldsymbol{x})$ given by Eq. 43 with $t_{\star}(\boldsymbol{x})=t_{e}(\boldsymbol{x})$ and compute its escape time $t_{e}^{\prime}\equiv t_{e}(\boldsymbol{x}_{i}^{\prime})$. The assumption is valid when, on average, $t_{e}^{\prime}-t_{e}\ll t_{e}$ for large $t_{e}$. We did this experiment in the following systems: open tent map with $a=3$ and $b=5$, standard map with $K=6$, and coupled Hénon maps with $D=2,4,6,8$ (see appendix VII.1) for $\Delta=1$. We observed that in all 6 cases, the average $t_{e}^{\prime}-t_{e}$ is smaller than 1 for all $t_{e}>2/\kappa$. These observations show that the assumption of Eq. 57 is valid for a broad class of chaotic systems.

A second assumption tested here is the self-similarity argument used in deriving Eq. 58. The self-similarity assumption we use in Eq. 58 is that the landscape is self-similar such that, irrespectively of the particular scale we choose (by decreasing $t_{\star}$), the properties of the zoomed phase-space remain the same. In particular, we are interested in checking that $t_{e}^{\prime}-t_{e}$ does not depend on $t_{e}$ when we choose a larger scale, i.e. when $t_{\star}$ is decreased by a constant. We thus repeat the experiment above but we decrease $t_{\star}$ to $t_{\star}(\boldsymbol{x})=t_{e}(\boldsymbol{x})-T$ with $T=0,1/\kappa,5/\kappa$ to check that increasing $T$ decreases the average $t_{e}^{\prime}-t_{e}$ without changing its independency with $t_{e}$. Our numerical tests in the same systems as before show that $t_{e}^{\prime}-t_{e}$ decreases with increasing $T$ and it remains independent of $t_{e}$, confirming our hypothesis.

The previous tests indicate that the proposal distribution should induce a constant acceptance rate when the Lyapunov exponent of the system, Eq. 74, is used. To confirm that this is the case, a flat-histogram simulation with an isotropic proposal distribution with width $\delta_{x}(\boldsymbol{x})=\delta_{x}(t_{e}(\boldsymbol{x}))$ given by Eq. 74 with $\lambda_{L}(t_{e})=\lambda_{L}$ was made. The results, Fig 10, reproduced from Ref. Leitão et al. (2013), confirm that proposing with $\lambda_{L}$ guarantees a constant acceptance, and that any other exponent in Eq. 74 fails to achieve so.

The above tests confirm that the derivation made in Sec. IV.3.2 holds for a paradigmatic strongly chaotic open system. These tests also present a major advantage of using the approach in this paper: it allows to test the assumptions made on each step, something that other approaches, such as the ones in Refs. Sweet et al. (2001); Kitajima and Iba (2011); Laffargue et al. (2013), do not explicitly allow.

The efficiency of the simulation is tested in Fig. 11 as a function of the maximal escape time considered, $t_{\max}$, $\tau(t_{\max})$ in the generic coupled Hénon maps defined by Eq. 69. It confirms the dramatic improvement of Metropolis-Hastings with the proposal derived in section II.2 over uniform sampling: the scaling is polynomial using importance sampling, and exponential in uniform sampling.

The derivation in Sec. IV.3.2, the tests presented above, and results in Figure 11, show the why and how importance sampling Metropolis-Hastings can efficiently sample long-living trajectories in strongly chaotic open systems. The proposal distribution should be applicable to strongly chaotic open systems more generally, as the approximations made are expected to be valid in other strongly chaotic systems.

VI Conclusions

VI.1 Summary of results

We have introduced a framework to sample rare trajectories in different classes of chaotic systems. It is based on the Metropolis-Hasting algorithm, a flexible and well-established Monte Carlo method suitable for the investigation of many numerical problems in chaotic dynamical systems (as shown in Sec. II). Our main contribution is a procedure (see Sec. IV) to construct the proposal step of the Metropolis-Hasting algorithm (which proposes a new state $\boldsymbol{x}^{\prime}$ given the current state $\boldsymbol{x}$) that ensures the efficiency of the sampling. The main arguments in the construction of this procedure are:

• (i)

set (in Sec. IV.1) as an heuristic goal to have a bounded (or constant) acceptance rate (23). This generalizes the traditional heuristic Viana Lopes (2006) used in Metropolis-Hastings, $E(\boldsymbol{x}^{\prime})-E(\boldsymbol{x})\sim 1$ and can therefore be used more generally in Metropolis-Hastings simulations.101010 For example, using this heuristics, in a Metropolis-Hastings flat-histogram in the Ising model (in this case $\boldsymbol{x}$ corresponds to a list of all spins, $E(\boldsymbol{x})$ corresponds to the energy of the configuration), when the configuration of the system is close to the maximum of the density of states (maximal energy, no magnetization), there is no need to flip just one spin at the time: one can flip all spins at once because, there, correlating such states brings no advantage to increase the acceptance rate (it is going to be accepted anyway), but it increases $r$.

• (ii)

introduce (in Sec. IV.2) an auxiliary quantity, the correlation time $t_{\star}$, that quantifies the similarity between any two states. It is motivated by the notion that the observables considered in section II.3 are computed over trajectories, and $t_{\star}(\boldsymbol{x},\boldsymbol{x}^{\prime})$ quantifies the time in which trajectories remain close to each other. We then showed how two proposal distributions, shift and neighborhood, can be used to control the average of $t_{\star}(\boldsymbol{x},\boldsymbol{x}^{\prime})$ over $\boldsymbol{x}^{\prime}$, $t_{\star}(\boldsymbol{x})$.

• (iii)

derive (in Sec. IV.2.1) an expression for the values of $t_{\star}(\boldsymbol{x})$ – Eq. 54 and 61 – which should be used in order to guarantees a constant acceptance. This is done for two observables of interest in chaotic systems (see Sec. II) – the escape time $E(\boldsymbol{x})=t_{e}(\boldsymbol{x})$ and the finite-time Lyapunov exponent (FTLE) $E(\boldsymbol{x})=\lambda_{t_{o}}(\boldsymbol{x})$ – and two target distributions of the Metropolis Hasting method (see Sec. III.1) – canonical and flat-histogram. These results are summarized in Tab. 1.

A successful application of our framework leads to an algorithm in which the number of samples to obtain an independent rare sample scales polynomially with the difficulty of the problem, as opposed to the exponential increase observed in traditional uniform sampling (see Fig. 9 and 11).

VI.2 Comparison to previous results

The importance of the proposal distribution has long been emphasized for Monte Carlo methods Roberts et al. (1997); Robert and Casella (2005) and for sampling chaotic systems Sweet et al. (2001); Dellago et al. (2002); Bolhuis et al. (2002); Tailleur and Kurchan (2007). The questions that remain open from these works, and that we tackle in our paper, are how the efficiency of the sampling method is related to the proposal distribution and how an efficient proposal can be constructed from assumptions about the system. The proposals derived in this paper can be mapped, under appropriate simplifications, to known results from the literature, specifically to the proposals derived Refs. Leitão et al. (2013, 2014) and the stagger part of the algorithm presented in Ref. Sweet et al. (2001). We make this connection explicitly in appendix VII.3. Another example of an algorithm that can be directly analyzed by the framework developed here is the precision shooting proposed in Ref. Grünwald et al. (2008); Grünwald (2009) and used in numerous applications of transition path sampling to chemical reactions Grünwald (2009); Rowley and Woo (2009); Eidelson and Peters (2012); Schaefer et al. (2014). The precision shooting method proposes a state $\boldsymbol{x}^{\prime}$ isotropically distanced from $\boldsymbol{x}$ by $\delta_{x}(\boldsymbol{x})$ given by Eq. 43, with $t_{\star}(\boldsymbol{x})=t_{o}$ where $t_{o}$ is the length of the trajectory. From our results we conclude that this proposal is sub-optimal because it over-correlates the proposed state. Another application of the results of section IV is in the algorithm Lyapunov Weighted Dynamics of Refs. Tailleur and Kurchan (2007); Laffargue et al. (2013), which uses a population Monte Carlo algorithm. As mentioned in these references (see also Ref. Wouters and Bouchet (2016)), there is a parameter $\varepsilon$ that controls how far new clones $\boldsymbol{x}^{\prime}$ should be distanced from the existing clone $\boldsymbol{x}$, and that it should be neither too small nor too large. This plays a role similar to $\delta_{x}(\boldsymbol{x})$ in section IV and an optimal $\varepsilon$ should therefore be related to our results.

In comparison to previous sampling methods in dynamical systems, including those mentioned above and others (e.g., ref. Sweet et al. (2001)), the distinguishing feature of our results is that they provide an explicit connection between the proposal distribution and the acceptance rate. This connection, which is typically absent in Monte Carlo methods more generally, is extremely powerful because failures of the algorithm can be related to violations of the hypothesis (about the method and dynamical system) that we used in our derivations. Such violations should then be understood, and this understanding can then be inserted back in this methodology to generate new methods adapted for that situation. We hope this process will increase the range of applicability of our framework to other classes of dynamical systems (e.g., non-hyperbolic systems Zaslavsky (2002); Altmann et al. (2013)) and observables $E$. Possible improvements of our results can be obtained considering anistoropic search domains, an idea that has shown to be essential in the case of finding chaotic saddles in systems with more than one positive Lyapunov exponents Sala et al. (2016).

VI.3 Discussion

The proposal distribution is a way of moving in the phase-space stochastically and how to select a new state $\boldsymbol{x}^{\prime}$ from a given state $\boldsymbol{x}$ is a general problem in different numerical techniques. Some of the most successful numerical algorithms in the literature, such as the golden section search or gradient descent, are essentially generic and efficient ways of selecting a new state. The results in Figs. 7, and  9 show that the proposal distribution strongly influences the computational cost of the different procedures, often irrespectively of the particular sampling procedure (canonical or flat-histogram) and problem (sampling or finding). This reinforces the notion that the proposal of the new tentative state from the current state is a crucial factor when developing numerical techniques for optimization and numerical integration. There we can expect that the main insights of this papers, e.g. the direct connection between fundamental properties of chaotic systems and the optimal proposal distribution, to be useful also for other problems (e.g., to the optimization problem in which one is interested in maximizing or minimizing the observable).

The development of a numerical algorithm requires compromising between how fast it solves a particular problem, and how it is able to solve different problems. One interesting aspect of the algorithms (i.e. proposal distributions) introduced in this paper is that even though they can be made very specific (e.g. propose with the FTLE of the trajectory, Eq. 43), they can also be made more general (e.g. propose using the Lyapunov exponent of the system, or the power-law proposal distribution, that does not use any specific information about the system). That is, more specificity requires more information (the FTLE of the trajectory) and makes the algorithm more efficient, and less information (only the Lyapunov of the system) makes the algorithm less specific, but also less efficient. This demonstrated capability of this methodology shows how it is not only useful to study a particular system on which some information about it is known, but also useful to situations on which less is known about the system. This does not mean that the methods apply to all problems, as there are important classes of chaotic systems on which some of the assumptions used in section IV are violated. For example, in non-hyperbolic systems Zaslavsky (2002); Altmann et al. (2013) trajectories may remain correlated for a long time, which implies that one cannot assume that after $t_{\star}$ the trajectories are independent. Nevertheless, because the framework was outlined in the form of adding known information about the system, it is possible that improved insights about a class of chaotic systems can be translated to a faster algorithm.

One advantage of the methodology presented here is that it is not restricted to specific observables $E(\boldsymbol{x})$. In principle, it can been used to construct proposal distributions to sample rare states in different observables $E$, $E(\boldsymbol{x})=t\lambda_{t}$, and $E(\boldsymbol{x})=t_{e}(\boldsymbol{x})$. While it remains unclear what is the precise class of observables for which our methodology allows to construct an efficient proposal distribution, the different derivations of the proposal distribution do provide insights on the observables for which a proposal distribution could be constructed from properties of the system. As argued at the end of Sec IV.3.1, observables computed as an average along the trajectories are similar to the FTLE and therefore the proposal distribution, derived in Eq. 56 should lead to efficient algorithms in these cases.

Altogether, our results reveal a fascinating interplay between the chaotic nature of some non-linear systems and the numerical techniques available to study rare events. It reinforces the idea that an efficient proposal requires information about the system (e.g., our derivation of the proposal distribution used the fact that “trajectories diverge exponentially” and that “the escape time function is a fractal-like function”). The analysis of this interplay allows to both better understand these systems and better understand these numerical techniques. This understanding opens perspectives to develop better techniques to numerical study rare trajectories and extreme events in non-linear systems more generally.

VII Appendices

VII.1 Appendix 1: Dynamical Systems

In this appendix we describe the different dynamical systems that we use throughout the paper.

VII.1.1 Skewed Tent Map

A paradigmatic example of a strongly chaotic system is the tent map Ott (1993), defined on $\Omega=[0,1]$ by

[TABLE]

where $a>1$ is a constant and $b\equiv a/(a-1)$ This map contains the main features of a chaotic system: it has a positive Lyapunov exponent ($\lambda_{L}=\frac{a\log(b)+b\log(a)}{a+b}$) and a positive measure. The finite-time Lyapunov exponent is given by

[TABLE]

where $i(x)$ is the number of times $x_{t}\in[0,1/a]$. Its distribution of the finite-time Lyapunov exponent for $U(\boldsymbol{x})=1$, $P(E)=G(E)$, can be computed analytically and is a binomial,

[TABLE]

VII.1.2 Standard map

The standard map considered here is defined by $\boldsymbol{x}\equiv(p,\theta)\in\Omega=[0,1]\times[0,1]$ that evolves in time according to

[TABLE]

We focus on the parameter $K=6$ which leads to a phase-space with no visible KAM islands. We also consider the leaked (open) version this map by introducing a hole into the system at $\Lambda=[0.1,0.1]$.

VII.1.3 Skewed Open Tent Map

The paradigmatic example of a strongly chaotic open system is the open tent map Lai and Tél (2010), defined on $\Omega=[0,1]$ by

[TABLE]

where $a>1$ and $b>a/(a-1)$. The state exits the system when it leaves the unit interval, i.e. $\Lambda=\mathbb{R}-\Omega$. This map contains the main features of an open chaotic system: it has a positive Lyapunov exponent ($\lambda_{L}=\frac{a\log(b)+b\log(a)}{a+b}$), an exponential decay of the escape time distribution,

[TABLE]

with $\kappa=-\log(1/a+1/b)$, and a conditionally invariant measure that is fractal with a (non-integer) fractal dimension $D_{0}$ given implicitly by

[TABLE]

VII.1.4 Coupled Open Hénon Map

As a generic example of a high-dimensional strongly chaotic open system, we consider a set of $d$ coupled Hénon maps on a ring, defined by a state $\boldsymbol{x}=(x_{1},y_{1},...,x_{d/2},y_{d/2})\in\Omega=\mathbb{R}^{d}$ where each individual map $(x_{i},y_{i})$ evolves according to

[TABLE]

for $i=1,...,d/2$, $d/2+1\equiv 1$, and with parameters $k=0.4$, $B=0.3$, $A_{1}=3$ (if $d>1$), $A_{d/2}=5$, and $A_{i}=A_{1}+(A_{d/2}-A_{1})(i-1)/(d/2-1)$. This choice of parameters ensures that a chaotic map is obtained in the $d=2$ case and corresponds to the map studied in Ref. Sweet et al. (2001) for $d=4$. The constraining region is $\Gamma=[-4,4]^{d}$ because it covers the chaotic saddle of the system, and $\Lambda=\Omega-\Gamma$, i.e. the trajectory leaves the system if the absolute value of any of the coordinates is higher than 4. The escape function $E(\boldsymbol{x})=t_{e}(\boldsymbol{x})$ for $d=4$ is represented in figure 2.

VII.2 Appendix 2: Efficiency of the uniform proposal

Here we show that the choice of sampling distribution does not necessarily decrease the scaling of the variance of the estimator. Our goal is to compute the average acceptance rate for a given escape time $t_{e}$ in the canonical ensemble with a uniform proposal distribution, $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})=1/|\Gamma|$. The acceptance is given by $a(\boldsymbol{x}^{\prime}|\boldsymbol{x})=\min\{1,\exp(-\beta(t_{e}(\boldsymbol{x}^{\prime})-t_{e}(\boldsymbol{x}))\}$ where $\beta<0$ is used to reach higher $E(\boldsymbol{x})=t_{e}(\boldsymbol{x})$. The acceptance of a state $\boldsymbol{x}$ is given by

[TABLE]

Because $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})=1/|\Gamma|$ and $\pi$ only depends on $t_{e}(\boldsymbol{x})$, $a(\boldsymbol{x})$ does not depend on $\boldsymbol{x}$, only on $t_{e}$: $a(\boldsymbol{x})=a(t_{e}(\boldsymbol{x}))$. The average acceptance rate at a given $t_{e}$, $A(t_{e})\equiv\mathbb{E}\left[a(\boldsymbol{x})|t_{e}\right]$, is given by

[TABLE]

Because $a(\boldsymbol{x})$ only depends on $t_{e}$, it can be pulled out of the integral, and thus $A(t_{e})=a(t_{e})$. Taking into account that $P(t_{e}^{\prime})=\kappa\exp(-t_{e}^{\prime}\kappa)$, the acceptance rate can be computed analytically by integrating Eq. 70 and leads to

[TABLE]

This shows that the acceptance rate decays exponentially with increasing $t_{e}$ (recall that $\beta<0$). Because low acceptance implies that the random walk stays on the same state for a long time, this leads to an exponential increase of the autocorrelation time $T(E)$ and therefore an increase of the variance in Eq. 22. This same argument applies to a flat-histogram simulation, where $\pi(\boldsymbol{x})\propto\exp(\kappa t_{e}(\boldsymbol{x}))$.

VII.3 Appendix 3: Simplified proposals

This section presents approximations that can be used to simplify both the implementation time and the computational cost of the proposals derived in Sec. IV. These approximations often reduce the general proposal derived in Sec. IV to particular proposals already found in the literature, and therefore explains such proposals in this wider context.

VII.3.1 Propose with the Lyapunov exponent in open systems

Eq. 43 requires computing $\lambda_{t_{e}(\boldsymbol{x})}(\boldsymbol{x})$, even though the main interest is in $t_{e}(\boldsymbol{x})$. This calculation requires using a numerical algorithm or multiply a product of matrixes Cvitanovic et al. (2016), both of which have an associated computational cost. A simplification to this proposal is to approximate $\lambda_{t_{e}(\boldsymbol{x})}(\boldsymbol{x})$ by the maximum of the distribution of FTLE with finite-time $t_{e}$, $\lambda_{L}(t_{e})$,

[TABLE]

This approximation is valid as long as $\lambda_{t_{e}}(\boldsymbol{x})$ is not on the tails of the distribution of FTLE with finite-time $t_{e}$, $P(\lambda_{t_{e}})$. A sampling distribution that only depends on $t_{e}$, $\pi(\boldsymbol{x})=\pi(t_{e}(\boldsymbol{x}))$, guarantees that states with the same $t_{e}$ are equally sampled, $P(\boldsymbol{x}|t_{e})=U(\boldsymbol{x})$, and therefore $P(\lambda_{t_{e}}(\boldsymbol{x})|t_{e})=P(\lambda_{t_{e}})$. Furthermore, the approximation of using the maximum of the distribution holds because the tails of $P(\lambda_{t_{e}})$ decay exponentially with increasing $t_{e}$ (see sec. II.1). Under these approximations, Eq. 43 can be simplified to

[TABLE]

Furthermore, $\lambda_{L}(t_{e})$ converges to the Lyapunov exponent of the system, $\lambda_{L}$, with increasing $t_{e}$. Therefore, a further simplification is to use the Lyapunov exponent of the system instead of $\lambda_{L}(t_{e})$ in Eq. 74,

[TABLE]

Furthermore, the $t_{\star}(\boldsymbol{x})$ we derived in Eq. 61, in the flat-histogram ensemble and with $\log G(t_{e})\propto-\kappa t_{e}$, can written as

[TABLE]

Defining the constant $\delta_{0}\equiv\Delta e^{-a\lambda_{L}/\kappa}$, we can write

[TABLE]

This equation is exactly the proposal derived in Ref. Leitão et al. (2013), and shows that the proposal derived in Sec. IV.3.2 generalises the proposal in Ref. Leitão et al. (2013).

VII.3.2 Adaptively estimate the Lyapunov exponent

Using the proposal distribution with $\lambda_{L}$ requires a priori knowledge of it, which typically is not available. This difficulty resembles the same problem that flat-histogram simulations have: $G(E)$ is required, but it is typically unknown a priori. This analogy motivates a Monte Carlo procedure that on the fly computes $\delta_{x}(t)$ that scales with $\lambda_{L}$.

Consider an hypothetical simulation with an isotropic proposal distribution (Eq. 41) with

[TABLE]

where $\sigma(t)$ is initially set to be $\sigma(t)=1$ for every $t$. Consider also that the simulation reached a state $\boldsymbol{x}$ with a high escape time (e.g. $t_{e}=t_{e}(\boldsymbol{x})=10/\kappa$). A proposed state, $\boldsymbol{x}^{\prime}=\boldsymbol{x}+\hat{\boldsymbol{h}}\sigma(t_{e})$, will most likely have a much lower escape time (e.g. $t_{e}(\boldsymbol{x}^{\prime})=1/\kappa$). From Eq. 43 and Eq. 61, this indicates that $\sigma(t_{e})$ is much higher than the “correct” proposal, $\delta_{x}(\boldsymbol{x})$, and therefore it should be reduced in the next proposal. The opposite is also true: when $\sigma(t_{e})$ is much smaller than $\delta_{x}(\boldsymbol{x})$, $t_{e}(\boldsymbol{x}^{\prime})=t_{e}(\boldsymbol{x})$ and it should be increased. This hypothetical simulation suggests that, in the same spirit as the Wang-Landau algorithm to approximate the density $P(t_{e})$, there is the possibility to approximate $\delta_{x}(\boldsymbol{x})$ using an update scheme that can be inserted in the Metropolis-Hastings algorithm, and that is given by the same algorithm as the Wang-Landau (see sec. III.3), but instead of updating $P_{WL}(t)$, it updates also $\sigma(t)$ Leitão et al. (2013):

[TABLE]

This update scheme generalises the Wang-Landau procedure to the proposal distribution. It is expected to converge to a function $\sigma(t)$ that decays exponentially with the Lyapunov exponent of the system, and a proposal distribution with a constant acceptance rate in a flat-histogram simulation. It was extensively tested in different systems (tent map, full chaotic standard map with a leak, Coupled Hénon map with different $D$s), see Refs. Leitão et al. (2013); Sala et al. (2016); Leitão (2016).

VII.3.3 Power-law proposal distribution

The proposal distributions derived in the previous sections requires some knowledge about the state and the system: $\lambda_{t_{o}}(\boldsymbol{x})$, $t_{e}(\boldsymbol{x})$ (in open systems), $\lambda_{L}$ of the system, and, in some situations, $G(E(\boldsymbol{x}))$. Another alternative to avoid computing $\delta_{x}(\boldsymbol{x})$ is to consider a proposal on which the time $t_{\star}$ that the two trajectories remain together is not imposed by a given $t_{\star}(\boldsymbol{x})$, but that is a uniformly random variable between $[0,t_{o}]$ (FTLE) or $[0,t_{e}(\boldsymbol{x})]$ (open systems) that is generated on each proposal. Some values of $t_{\star}$ will be far from the optimal $t_{\star}(\boldsymbol{x})$ and the corresponding $\boldsymbol{x}^{\prime}$ will be rejected or it will be too close from $\boldsymbol{x}$, but others $t_{\star}$ will still be close from the optimal $t_{\star}(\boldsymbol{x})$ and therefore useful.

Having a uniformly distributed correlation $t_{\star}$ still requires computing $\delta_{x}(\boldsymbol{x})$ in Eq. 43, which requires $\lambda_{t_{o}}(\boldsymbol{x})$. In the case $\lambda_{t}(\boldsymbol{x})$ is unknown (e.g. in open systems one could be only interested in the escape time and therefore not compute $\lambda_{t}(\boldsymbol{x})$), one may further approximate it by an uniform distribution between two extremes. Because the product of two uniformly distributed random variables is also uniformly distributed, this leads to a proposal distribution where $\delta_{x}(\boldsymbol{x})$ is given by $\exp(-U(a,b))$ where $a$ and $b$ are free parameters. By standard transformation of variables, this leads to a scale $\delta_{x}(\boldsymbol{x})$ that is power-law distributed and given by

[TABLE]

where $\delta_{\min}$ and $\delta_{\max}\approx|\Gamma|$ are two free parameters. The $\delta$ in $\boldsymbol{x}^{\prime}=\boldsymbol{x}+\boldsymbol{h}\delta$ is an half-normal distribution with a scale $\delta_{x}$, but since this scale is now power-law distributed, it is no longer necessary to use the half-normal distribution altogether; instead, it is possible to just use the power-law proposal distribution where $\delta$ is drawn from $P(\delta_{x})$ in Eq. 80, i.e. $|\boldsymbol{x}^{\prime}-\boldsymbol{x}|$ is power-law distributed according to Eq. 80. Without the half-normal distribution the proposal distribution no longer depends on $\boldsymbol{x}$ and therefore $g(\boldsymbol{x}^{\prime}|\boldsymbol{x})/g(\boldsymbol{x}|\boldsymbol{x}^{\prime})=1$.

The stagger part of the algorithm of Ref. Sweet et al. (2001) proposes exactly with a scale given by Eq. 80 and the argumentation above explains why the proposal distribution used in Ref. Sweet et al. (2001) to find states with high-escape time $t_{e}$ is reported to work well: it is a proposal distribution that proposes $\boldsymbol{x}^{\prime}$ correlated with $\boldsymbol{x}$ with a correlation $t_{\star}$ that is uniformly distributed, which eventually proposes $\boldsymbol{x}^{\prime}$ with the optimal correlation $t_{\star}(\boldsymbol{x})$. To confirm this explanation, let us consider a flat-histogram simulation with a power-law proposal distribution on the open tent map and consider the measurement of $\mathbb{E}\left[\log\delta_{x}|A^{*},t_{e}\right]$, where $A^{*}$ is the condition $\varepsilon<a(\boldsymbol{x}^{\prime}|\boldsymbol{x})<1-\varepsilon$ (of bounded acceptance). Under the above argumentation, the scale $\log\delta_{x}$ that contributes to a bounded acceptance is given by $-\lambda_{t_{e}}t_{e}$, per Eq. 77. This is confirmed by numerical simulation, shown in figure 14, and was obtained also for the problem of finding rare states, as reported in Ref. Sala et al. (2016). This result, combined with the derivation of $t_{\star}(\boldsymbol{x})$, explains the success of the proposal (the stagger) used in Ref. Sweet et al. (2001) from basic notions of chaotic systems and numerical methods.

VII.4 Appendix 4: Algorithmic description

It is useful to summarise the different proposals in an algorithmic form so they can be easily implemented (see Ref. Leitão (2017) for our codes). In all cases, the proposal requires the current state of the random walk, denoted by $\boldsymbol{x}$.

VII.4.1 FTLE in closed systems

Generate a unitary vector $\hat{\boldsymbol{\delta}}$ in $D$ dimensions 2. 2.

Compute $t_{\star}=t_{o}-\frac{\log(a)}{d\log\pi(E)/dE(E)}\frac{1}{|\lambda_{L}-\lambda_{t_{o}}(\boldsymbol{x})|}$, Eq. 54. 3. 3.

Compute $\delta_{x}=\Delta\exp(-t_{\star}\lambda_{t_{o}}(\boldsymbol{x}))$, Eq. 43 4. 4.

Generate a random number $\delta$ from a normal distribution with mean 0 and variance $\delta_{x}^{2}$ 5. 5.

Make $\boldsymbol{x}^{\prime}=\boldsymbol{x}+\hat{\boldsymbol{\delta}}|\delta|$

where $\Delta\in\mathbb{R}$ is free parameter (e.g. $\Delta=0.1$) and $a$ is the chosen average acceptance (e.g. $0.5$). For example, in a canonic ensemble with parameter $\beta$, $d\log\pi(E)/dE(E)=\beta$ and therefore $t_{\star}=t_{e}(\boldsymbol{x})-\frac{\log(a)}{\beta}\frac{1}{|\lambda_{L}-\lambda_{t_{o}}(\boldsymbol{x})|}$. The value of $\lambda_{L}$ can be estimated using e.g. the first samples of the random walk. In the flat-histogram ensemble, $\lambda_{L}$ is the maximum of $G(\lambda_{t_{o}})$ and $\log\pi(E)/dE(E)=\log G(E)/dE(E)$, or by an approximation of it, e.g. using $G_{WL}$ of the Wang-Landau algorithm.

VII.4.2 Open systems

Generate a unitary vector $\hat{\boldsymbol{\delta}}$ in $D$ dimensions 2. 2.

Compute $t_{\star}=t_{e}(\boldsymbol{x})-1/\kappa-\frac{a-1}{d\log\pi/dt_{e}}$, Eq. 61 3. 3.

Compute $\delta_{x}=\Delta\exp(-t_{\star}\lambda_{t_{o}}(\boldsymbol{x}))$, Eq. 43 4. 4.

Generate a random number $\delta$ from a normal distribution with mean 0 and variance $\delta_{x}^{2}$ 5. 5.

Make $\boldsymbol{x}^{\prime}=\boldsymbol{x}+\hat{\boldsymbol{\delta}}|\delta|$

where $\delta_{0}$ is a free parameter (e.g. $\delta_{0}=0.1$). Both $\lambda_{t_{e}}(\boldsymbol{x})$ and $t_{e}(\boldsymbol{x})$ are required by the proposal and both can be computed during the same evolution of the system: $t_{e}(\boldsymbol{x})$ is the time until the trajectory enters the exit region $\Lambda$, $\lambda_{t_{e}}(\boldsymbol{x})$ is the FTLE of this trajectory.

Acknowledgements.

JCL was funded by FCT (Portugal), Grant NO. SFRH/BD/90050/2012. EGA thanks J. Wouters and T. Tél for helpful discussions.