Taming chaos to sample rare events: the effect of weak chaos

TL;DR

This paper explores strategies to efficiently sample rare events in weakly chaotic dynamical systems by controlling non-hyperbolicities, extending previous methods used for strongly chaotic systems.

Contribution

It introduces new approaches to overcome limitations posed by non-hyperbolicities in weakly chaotic systems within a Monte Carlo framework.

Findings

Identified how non-hyperbolicities hinder sampling efficiency.

Proposed strategies to improve sampling in weakly chaotic systems.

Validated methods on low-dimensional chaotic maps.

Abstract

Rare events in non-linear dynamical systems are difficult to sample because of the sensitivity to perturbations of initial conditions and of complex landscapes in phase space. Here we discuss strategies to control these difficulties and succeed in obtainining an efficient sampling within a Metropolis-Hastings Monte Carlo framework. After reviewing previous successes in the case of strongly chaotic systems, we discuss the case of weakly chaotic systems. We show how different types of non-hyperbolicities limit the efficiency of previously designed sampling methods and we discuss strategies how to account for them. We focus on paradigmatic low-dimensional chaotic systems such as the logistic map, the Pomeau-Maneville map, and area-preserving maps with mixed phase space.

Equations

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

$$\pi (x)=\pi (E_x)\propto e^{-\beta E_x},$$

$$\pi (x)\propto \frac{1}{P(E_x)}\text{ for }{E}_x\in [E_{min},E_{max}],$$

$$g(x|x^{\prime })\approx g(x^{\prime }|x),$$

$$\mathbb{E}\left[\frac{\pi (E_{x^{\prime }})}{\pi (E_x)}|x\right]={\int }_{\Gamma }\frac{\pi (E_{x^{\prime }})}{\pi (E_x)}g(x^{\prime }|x)d{x}^{\prime }=a,$$

$$\frac{\pi (E_{x^{\prime }})}{\pi (E_x)}=1+\frac{d\log \pi (E_x)}{dE}(E_{x^{\prime }}-E_x).$$

$$\mathbb{E}\left[E_{x^{\prime }}-E_x|x\right]=\frac{a-1}{d\log \pi (E_x)/dE}.$$

$${\delta }_x(x)=\Delta e^{-{\lambda }_{t_{\star }}(x)t_{\star }(x)},$$

$$x_{t+1}=F(x_t)=4x_t(1-x_t).$$

$$\left|x_{t_i+1}^{\prime }-x_{t_i+1}\right|\approx \left|\frac{dF}{dx}(x_{t_i})\right|{\delta }_0\approx 0.$$

$$x_{t_i+1}^{\prime }-x_{t_i+1}\approx \frac{1}{2}\left|\frac{d^{2}F}{d{x}^2}(x_{t_i})\right|{\delta }_0^2=4{\delta }_0^2.$$

$$x_{t+1}=F(x_t)=x_t+x_t^2\mathrm{mod}1.$$

$$\ell =\frac{1}{2}\left(-1+\sqrt{3+2\sqrt{5}}\right),$$

$${\delta }_x(x,t_{\star })=\Delta e^{-{\lambda }_{t_{\star }}{t}_{\star }}.$$

$$\log t_e(x^{\prime })⪆\log t_e(x)+\log q$$

$$\Delta E(x)=\log t_e({\epsilon }^{\prime })-\log t_e(\epsilon )=\log \epsilon -\log {\epsilon }^{\prime }.$$

$${\delta }_i(x)=\left(e^{-\Delta E_x}-1\right)\frac{1}{t_e(x)}.$$

$${\delta }_x(x)=\Delta e^{-{\lambda }_{t_i}(x)t_i}{\delta }_i(x)=\Delta e^{-{\lambda }_{t_i}(x)t_i}\left(e^{-\Delta E_x}-1\right)\frac{1}{t_e(x)}.$$

$${\delta }_x(x)=\Delta \left(e^{-\Delta E_x}-1\right)e^{-{\lambda }_{t_e}(x)t_e}{t}_e(x),$$

$$F(p,\theta )=\{\begin{array}{l}p+K/(2\pi )\sin (2\pi \theta )\mathrm{mod}1\\ \theta +p+K/(2\pi )\sin (2\pi \theta )\mathrm{mod}1\end{array}.$$

$$\frac{g(x|x^{\prime })}{g(x^{\prime }|x)}=\frac{{\delta }_x(x)}{{\delta }_x(x^{\prime })}\exp \left[-\frac{\pi |x^{\prime }-x{|}^2}{4{\delta }_x(x{)}^2}\left(1-\frac{{\delta }_x(x{)}^2}{{\delta }_x(x^{\prime }{)}^2}\right)\right]$$

$$f(R)=R\exp \left[-\left(1-R^2\right)\pi /4\right].$$

$$P(\vec{\delta }|x)=P(\hat{\delta }|x)P(\delta |x),$$

$$P(\delta |x)=\frac{\sqrt{2}}{\sqrt{\pi {\delta }_x(x{)}^2}}{e}^{-\frac{\pi {\delta }^2}{4{\delta }_x(x{)}^2}}\text{ for }\delta >0.$$

Full text

**Forecast in chaotic dynamical systems requires the evolution of an ensemble of trajectories, which can quickly lead to very different outcomes. If the choice of the ensemble is compatible with our knowledge of the current state of the system, we can associate the probability of an event with the fraction of initial conditions for which it occurs. Computationally, this strategy is efficient to determine the most likely events but it struggles to compute rare events, e.g., those at the tail of the distribution of an observable of interest. The importance of such extreme events is that they often cause the largest impact and, due to the chaoticity of the system, they can not be easily anticipated. The development of efficient sampling methods, as aimed in this paper, is crucial in these cases because they are able not only to find trajectories leading to extreme events but also to estimate their probability (in the original ensemble). **

I Introduction

We are interested in sampling rare events in chaotic dynamical systems, a problem that has been the subject of different approaches in the recent years dellago2002transition ; tailleur2007probing ; geiger2010identifying ; iba2014multicanonical ; wouters2016rare . Given an initial condition $x$ in a $d$-dimensional phase space $\Omega$, $x\in\Omega\subset\mathbb{R}^{d}$, the dynamical system $F$ evolves it $x(t)=F^{t}(x(0))$ until (at time $t_{o}$) an observable $E_{x}=E(x(t_{o}))\in\mathbb{R}$ is measured. A rare event corresponds to an observable $E$ at the tail of the distribution $P(E)$ obtained from an ensemble of initial conditions selected in $\Gamma\subset\Omega$ according to a probability measure $\mu$ such that $\int_{\Gamma}d\mu=1$, e.g., $\mu$ can be simply the phase space volume (uniform distribution) or the natural measure of the dynamical system. Since the dynamics is deterministic, the variability of the events $E$ and the need for its statistical description are solely due to the sensitive to variations in the initial conditions and not due to an intrinsic random dynamics.

The importance-sampling method bucklew2013introduction we construct in this paper samples initial conditions $x\in\Gamma$ with probabilities different from $\mu$ in order to obtain more samples at the tail of $P(E)$. However, our goal is still to be able estimate the probability of the event $P(E)$ of our original problem (original ensemble of initial conditions). The strict determinism of the (chaotic) dynamics is crucial for the design of efficient sampling methods. On the one hand, the lack of intrinsic randomness poses difficulties to traditional methods, e.g., it is not possible to distinguish between trajectories based on different noise realizations and therefore cloning of trajectories requires more refined procedures tailleur2007probing ; wouters2016rare . On the other hand, determinism can be explored in order to efficiently search for trajectories leading to rare events.

The goal of this manuscript is to show how to construct efficient Markov Chain Monte Carlo methods to sample rare trajectories of chaotic dynamical systems. The key ingredient is to use information about the chaoticity of the last sampled trajectory to construct a proposal distribution that efficiently finds a new trajectory of interest. First we review our previously proposed  (leitao2017importance, ) approach (in Sec. II), which has been successful in different problems involving strongly chaotic systems leitao2017importance ; leitao2013monte ; gupta ; tapias . We then focus (in Sec. III) on deviations from strong chaos and how they pose challenges for the application of the previously developed methods.

II Chaos and Metropolis-Hastings methods

Sampling methods typically exploit the fact that if a given trajectory of interest is found – $x$ with a rare $E_{x}$ – this can be used to find other trajectories of interest – $x^{\prime}$ with $E_{x^{\prime}}\approx E_{x}$. The essential step for the success of the method is to be able to choose the right proposal distribution of $x^{\prime}$ given $x$ – denoted as $g(x^{\prime}|x)$ – that guarantees that $x^{\prime}$ will likely lead to a “good” $E_{x^{\prime}}$. The most natural choice is to correlate trajectories $x$ and $x^{\prime}$ to obtain $E_{x^{\prime}}$ sufficiently close to $E_{x}$. This is an heuristics often used in Statistical Physics, implicit in the choice of minimal/local proposals, e.g., in spin systems (single spin flip) Lopes.phd2006 , in random networks (single link exchange) Fischer2015 , and in proteins Grassberger1997 . Since our phase space is continuous it is not clear what a minimal local proposal would be in our case. In this sense, the challenge we address in this paper is to formalize what sufficiently close means and to construct an efficient proposal $g(x^{\prime}|x)$ that achieves it.

The general ideas sketched above are valid for broad classes of sampling methods, but here we focus on a Metropolis-Hastings setting NewmanBarkemaBook ; RobertCasellaBook , in line with our previous works revised in leitao2017importance . In this setting, starting from $x\in\Gamma$ a new state $x^{\prime}\in\Gamma$ is proposed according to $g(x^{\prime}|x)$. This move can be accepted – the new trajectory $x^{\prime}$ is sampled and the procedure is repeated from $x^{\prime}$ – or rejected – the trajectory $x$ is sampled again and the procedure is repeated from $x$ with an independent sample from $g(x^{\prime}|x)$. This procedure is repeated $n$ times, leading to $n$ (correlated) samples. If the proposal is ergodic – all $x\in\Gamma$ have a non-zero probability to be sampled for $n\rightarrow\infty$ – and the acceptance given by

[TABLE]

the sampled trajectories $x^{\prime}$ will approach $\pi(x)$ for $n\rightarrow\infty$ RobertCasellaBook . The sampling distribution $\pi(x)$ can be chosen at will. A popular choice is the canonical distribution NewmanBarkemaBook

[TABLE]

where different regions of the distribution $P(E)$ are sampled when the parameter $\beta$ is varied. Another popular choice is the multi-canonical (flat-histogram) distributionBerg1991

[TABLE]

which can be computed (in case $P(E)$ is unknown) also through the Wang-Landau method Wang2001 .

The crucial step to implement a Metropolis-Hastings algorithm to sample chaotic trajectories is the construction of an efficient proposal $g(x^{\prime}|x)$ distribution. The following three steps can be used to achieve this leitao2017importance :

The goal is to bound the acceptance (1) by making $E_{x}^{\prime}$ and $E_{x}$ sufficiently close to each other. Assuming

[TABLE]

the key remaining term in the acceptance is the ratio $\pi(E_{x^{\prime}})/\pi(E_{x})$ . We can thus express our heuristic more formally by fixing the expectation of this ratio over all possible $x^{\prime}$

[TABLE]

where $0<a\leq 1$ is a constant (ideally, the constant acceptance rate). Since the proposal achieves a small variation of $E$, $\pi(E_{x^{\prime}})$ can be expanded in Taylor series around $E_{x^{\prime}}=E_{x}$ as

[TABLE]

Introducing Eq. (6) in Eq. (5) we obtain an explicit condition

[TABLE]

• 2.

The next step is to compute the correlation time $t_{\star}$ needed for the trajectories to be close to each other in order to achieve condition (7). We assume that the observables $E$ of interest are built throughout the $t_{o}$ times step of the trajectory so that trajectories that remain close (correlated) in the phase space lead to similar observables $E$. The correlation time $t_{\star}$, $0\leq t_{\star}\leq t_{o}$ is the time the two trajectories remain “close” to each other, i.e., within a distance $\Delta$ that is smaller than the expected distance between two randomly chosen trajectories in $\Omega$. We assume that, in practice, the two trajectories are identical until $t_{\star}$, i.e. $x(t)=x^{\prime}(t)$ for $0\leq t\leq t_{\star}$,and independent for $t>t_{\star}$, i.e. $x^{\prime}(t)$ is sampled according to $\mu$ for $t>t_{\star}$. Explicit expressions for $t_{\star}$ have been derived for the escape time and finite time Lyapunov exponent (see Ref. leitao2017importance ) and for the dispersion in spatially extended (diffusive) systems (see Ref. tapias ).

• 3.

Once the $t_{\star}$ that guarantees condition (7) is known for a given problem (i.e., for a given observable $E$ and distribution $\pi$) we can generate trajectories $x^{\prime}$ from one of the following two procedures: shifting the trajectory by a time $t_{\star}$ backward/forward using the dynamics $x^{\prime}=F_{\pm t_{\star}}(x)$ (shift proposal); or proposing $x^{\prime}$ on a neighborhood of size $\delta_{x}$ around $x$ (local proposal). The choice of $\delta_{x}$ in a chaotic system is such that the trajectories should remain close to each other up to a time $t_{\star}$ despite the exponential divergence of nearby trajectories, and thus

[TABLE]

where $\lambda_{t}(x)$ is the largest finite-time Lyapunov exponent (FTLE) of the trajectory (initiated in position $x)$ and $\Delta$ is a constant of the order of $|\Gamma|$. Depending on the problem, $\lambda_{t_{\star}}$ can be approximated by $\lambda_{t_{o}}$ or (more strongly) by $\lambda_{t\rightarrow\infty}$ (the largest Lyapunov exponent of the system).

The construction proposed above has been successfully implemented to obtain efficient Monte Carlo methods in different problems leitao2017importance , including calculations of finite-time Lyapunov exponents in N coupled oscillators (with N up to $1024$) gupta and the computation of trajectories with high dispersion in diffusive systems such as the Lorentz gas tapias . Here we focus on violations of the simplifying assumptions and approximations made above. Violations of strong chaos and uniform hyperbolicity are typical in chaotic dynamical systems, and our interest is to investigate how they affect our sampling methods and how our methods can be modified to account for them. We denote by weak chaos the chaotic dynamics observed in systems that violate the simplifying hypothesis of uniform exponential divergence of initial conditions used above, including systems with marginal stable points and Hamiltonian systems with mixed phase space OttBook ; Zaslavsky2002 .

III Effect of weak chaos on the sampling algorithm

In each of the three subsections below we consider a simple dynamical systems, with increasingly important (generic) weakly-chaotic features, that violate some of the simplifying assumptions used above.

III.1 Logistic map

A crucial step in our derivation above is that the width of the local proposal, $\delta_{x}(x)$ in Eq. (8), guarantees that the trajectory starting from $x^{\prime}$ is within $\Delta$ of the trajectory starting from $x$ up to time $t_{\star}$. This approximation was based on the assumption of exponential divergence of nearby trajectories, using the maximum FTLE $\lambda_{t}(x)$. The example we consider below shows how this assumption can be violated due to the coexistence of regions with positive and negative FTLE. Consider the logistic map, defined on $\Omega=[0,1]$ by

[TABLE]

Any irrational initial condition $x$ ergodically fills $\Omega$ following $d\mu=\frac{1}{\pi\sqrt{x(1-x)}}dx$. The FTLE in this simple system can fluctuate considerably. The Lyapunov exponent is positive, but finite time estimations can be negative because $|dF/dx|=|4-8x|<1$ for $3/8<x<5/8$. Particularly problematic are trajectories that come close to $x=1/2$, where $|dF/dx|=0$. As a consequence, the distribution of FTLE of this map, $P(\lambda_{t})$, has negative values for any finite $t$ Prasad1999 , which implies that there is a non-zero measured set where $\lambda_{t}(x)<0$ for all $x$. For any $x$ in this set, and any fixed $\Delta$ in Eq. (8), increasing $t_{\star}$ leads to a proposal with a width larger than 1 (the phase space size). This implies that, for large $t_{\star}$, $x^{\prime}$ is approximately drawn uniformly from $[0,1]$. In this situation, it is never expected that the two trajectories $x^{\prime}$ and $x$ are close within $\Delta$ up to time $t_{\star}$, violating our initial assumption.

The crucial violation here is related to the fact that the Lyapunov exponent corresponds to a linear (first-order) term that correctly describes the divergence of two trajectories in time in the limit that their initial separation goes to 0. The growth of the divergence of two trajectories separated by a finite distance $\delta_{x}(x)$ is not necessarily well described by $\exp(-\lambda_{t}(x)t)$ (see Ref. vulpiani2013 for more on finite-size Lyapunov exponents). This violation is more evident in trajectories $x$ where $\lambda_{t}(x)<0$. To see why, consider a trajectory starting at $x$ on which at time $t_{i}$, $x_{t_{i}}\approx 1/2$, and consider that all other $x_{t}$ are not close to $1/2$. Since the map is chaotic, up to $t_{i}$, trajectories starting from $x^{\prime}$ distanced from $x$ by $\delta_{0}\exp(-\lambda_{t_{i}}(x)t_{i})$ were approximately, at time $t_{i}$, within $\delta_{0}$ of $x_{t_{i}}$. However, at that particular point $x_{t_{i}}$, the derivative is zero and thus the first order approximation predicts that the states $x_{t_{i}+1}$ and $x_{t_{i}+1}^{\prime}$ will be arbitrarily close to each other:

[TABLE]

However, this prediction is might not be accurate because when the first order term is zero, the second order term is non-zero and dominates:

[TABLE]

We confirmed numerically the appearance of the behaviour described above, which is known as glytch ChaosBook . We consider two trajectories initially separated by $\Delta\exp(-\lambda_{t}(x)t)$, and we compare the distance in time, $\delta_{t}\equiv|x_{t}^{\prime}-x_{t}|$, with the distance $\delta_{x}(t)=\Delta\exp(\lambda_{t}(x)t)$ expected based on the first order approximation (10). A violation of the assumption happens when the ratio $r(t)\equiv\delta_{t}/\delta_{x}(t)$ is different from $1$. Figure 1 shows a representative example of this simulation, which confirms that the ratio $r(t)$ can change abruptly, becoming orders of magnitude different from $1$ (indicating that the distance $|x_{t}^{\prime}-x_{t}|$ is much larger than the expected distance given by $\delta_{x}(t)$). This could be fixed by decreasing the initial distance $\Delta$, but the crucial point here is that $\Delta$ strongly depends on the particular $x$ and can be orders of magnitude different for different $x$ (e.g. one with $\lambda_{t}>0$ vs. one with $\lambda_{t}<0$). In other words, the assumption that is violated in the logistic map is that there is a $\Delta$ independent of $x$ that makes $\delta_{x}(x)$ in Eq. (8) to guarantee a correlation time $t_{\star}$ between any two trajectories. The results above do not imply that Metropolis-Hastings cannot be used in systems where $\lambda_{t}<0$ for some states, they imply that Eq. (8) has to be extended, e.g., to make $\Delta$ dependent on $x$. In the next section we consider a similar issue, arising when periodic orbits show zero Lyapunov exponents.

III.2 Pomeau-Maneville map

One important approximation in the derivation of $t_{\star}(x)$ in point 2. of Sec. II above is that, for $t>t_{\star}$, $x_{t_{\star}}^{\prime}$ is independent of $x_{t_{\star}}$ . This approximation was based on the notion that trajectories diverge exponentially and thus two trajectories are separated by $\Delta\approx 1$ at $t=t_{\star}$ will rapidly become independent of each other. This approximation is naturally violated when the (local) divergence of nearby trajectories is not exponential, and our goal here is to explore the consequences of this violation to our sampling method.

Let us analyze one simple one-dimensional system where non-exponential divergence is present, the Pomeau-Maneville map defined in $\Omega=\Gamma=[0,1]$ by

[TABLE]

This map is a model for intermittency, a phenomenon on which trajectories irregularly alternate between regular and chaotic motion OttBook . The intermittency in this system appears because $\frac{dF}{dx}=1+2x$ and thus the fixed point $x=0=F(0)$ is a non-hyperbolic point. The ergodic invariant measure $\mu$ is non-normalizable as it diverges at the fixed point as $d\mu\sim\frac{1}{x}dx$ Thaller1995 . More general Pomeau-Maneville maps consider a generic power $z$ instead of $2$ in Eq. (12) Niemann . From this point of view, the case treated here ($z=2$) is special because it lies at the border between normalizable and non-normalizable $\mu$ and it would be interesting to generalize our results to $z\neq 2$.

The observable we are interested in is the time a trajectory takes to escape an open Pomeau-Maneville map, achieved leaking Altmann2013 the map by adding an exit region $\Lambda$ so that trajectories $x\in\Lambda$ are removed ($t_{e}=0$). Choosing $\Lambda=[\ell,1]$ for the map (12), with

[TABLE]

ensures that the function relating the escape time $t_{e}$ to the initial condition $x$ can be described by a symbolic sequence with forbidden sequences, as shown in Fig. 2. The distribution of escape times $P(t_{e})$ is known to have a power-law tail $P(t_{e})\sim t_{e}^{-\alpha}$ with an exponent $\alpha=2$. This implies that the average escape time $\langle t_{e}\rangle$ diverges, a strong form of intermittency (or stickiness at the origin). Since the escape time varies over orders of magnitude, it is natural to consider as an observable the logarithm of the escape time $E_{x}=\log t_{e}(x)$. The distribution of $E_{x}=\log t_{e}(x)$ is then exponential with an exponent $\alpha^{\prime}=\alpha+1=-1$. Our interest is to estimate $P(t_{e})$ and sample trajectories at the tail of this distribution.

Qualitatively, a typical long living trajectory can be pictured by a trajectory that, for a time $t_{chaos}$, behaves as if it was a chaotic trajectory, and that at some time, denoted here as a time $t_{i}$, is injected close to the non-hyperbolic point $x=0$. The trajectory then spends a long time $t_{stick}$ close to [math], until it eventually leaves the region, returning to a chaotic movement. An example of such a trajectory is shown in Fig. 3a.

We now investigate how the recipe from Sec. II can be used to construct an efficient proposal distribution for the Pomeau-Maneville map, focusing on the main differences between this map and other strongly chaotic system for which the recipe has worked in the past. The two major differences here are: a) the observable is $E=\log t_{e}$, instead of $t_{e}$; and b) there is a non-hyperbolic fixed point at $x=0$. We obtain numerical insights on this problem by starting from a trajectory $x$ with a high escape time $t_{e}(x)$ and searching for different $x^{\prime}$, obtained adding to $x$ a small perturbation of typical size $\delta_{x}$ (see Appendix A for details) given by

[TABLE]

This equation is similar to Eq. (8), but here instead of using a theoretically derived $t_{\star}$ we use it as a free-parameter that defines a scale. Using this proposal, we measure $P(\log t_{e}^{\prime}-\log t_{e}|x)$ for different $t_{\star}$. The goal is to test the hypothesis that proposing with Eq. (8) guarantees that the two trajectories remain, on average, close together up to time $t_{\star}$. If the assumption holds in this system, choosing $t_{\star}=qt_{e}$, with $0<q<1$, would imply that on average $t_{e}(x^{\prime})\gtrapprox t_{\star}=qt_{e}(x)$. In terms of the logarithm, this would imply that

[TABLE]

The outcome of a numerical experiment that implements the ideas above is shown in Fig. 3. We focus on a trajectory $x$ with an escape time $t_{e}(x)=16458$, or $\log t_{e}\approx 10$ (Fig. 3, upper panel). The variation in the observable $P(\log t_{e}^{\prime}-\log t_{e}|x,t_{\star})$ for different values of $t_{\star}$ shows (Fig. 3, middle panel) that, independently of $t_{\star}$, most trajectories show a $\log t_{e}^{\prime}-\log t_{e}$ much smaller than the expected value from Eq. (15) (e.g., for $q=0.5$, $t_{\star}=0.5t_{e}$, we would expect $\log t_{e}(x^{\prime})-\log t_{e}(x)\gtrapprox-0.7$). In fact, independent of $t_{\star}$, almost 50% of all trajectories shows $\log t_{e}(x^{\prime})-\log t_{e}(x)\approx-6.5$. This shows that the assumption that Eq. (8) guarantees that the states are close up to $t_{\star}$ is violated here due to the non-hyperbolic nature of the point $x=0$. Still, our result does indicate a dependence of $\log t_{e}(x^{\prime})-\log t_{e}(x)$ on $t_{\star}$, which suggests that it may still be possible to derive a distance between $x^{\prime}$ and $x$ that leads to a bounded acceptance. To investigate this possibility, we repeat the approach done previously for the logistic map (in Fig. 1) and plot (in Fig. 3, lower panel) the expected divergence given by the first order term of the Taylor expansion with the actual distance between the trajectories, for different trajectories $x^{\prime}$ generated with $t_{\star}=0.9t_{e}$. We find that there are many trajectories that largely deviate from $x$ at the time $t_{i}$ when the trajectory is injected close to the critical point $x\approx 0$. Half of the trajectories quickly exit the system (being responsible for the high peak around $\log t_{e}(x^{\prime})-\log t_{e}(x)\approx-6.5$), they correspond to points that in Fig. 3 have a low escape time $t_{e}(x^{\prime})$.

The numerical observations reported above can be understood analytically by focusing on the injection of the chaotic trajectory into the trapping point $x=0$. This injection happens around the pre-image of $x=0$, which is $x=1/\phi$, where $\phi$ is the golden ratio. At $1/\phi$, the map $F$ in Eq. (12) is discontinuous. For a trajectory to be long living, it must approximate [math], which requires its pre-image to be very close but larger than $1/\phi$, i.e., $x=1/\phi+\epsilon$ (with $\epsilon>0$). The proposal distribution is a normal distribution around $x$, and therefore half of the trajectories can be at $x<1/\phi$. These will not be mapped close to [math] and therefore most likely they quickly leave the system. The orbits mapped $\varepsilon$ close to [math] evolve initially as $x_{t}\approx\varepsilon+t\varepsilon^{2}$ (in order $\varepsilon^{2}$). The escape time $t_{e}(x)$ is proportional to the time $x$ takes from $\varepsilon$ to $2\varepsilon$ and therefore $t_{e}\propto 1/\varepsilon$. The linear divergence close to $x\approx 0$ implies also that the total divergence of such trajectory is $D=\sum_{i=0}^{t_{e}(\varepsilon)}\log(1+2x_{i})\approx\sum_{i=0}^{t_{e}(\varepsilon)}2x_{i}\propto t_{e}^{2}$. We performed numerical simulations that confirm these two scaling and fix the pre-factors as $t_{e}=1/\varepsilon$ and $D=3.3t_{e}^{2}$.

We now use the results discussed above to obtain an expression for $\delta_{x}(x)$ that guarantees a given variation of $\log t_{e}(x^{\prime})-\log t_{e}(x)$ – as required by Eq. (7) – and that can thus be used to construct an efficient proposal distribution $g(x^{\prime}|x)$. Let $\varepsilon(x)$ be the smallest distance of $x_{t}$ from 0 for all $t=1,...,t_{e}(x)$, and let $t_{i}(x)+1$ be the time at which this happens (i.e., at $t_{i}(x)$, $x$ is injected close to $x=0$ ). Assuming that there are no re-injections, the time from $t_{i}$ until the state crosses $1/\phi$, $t_{e}^{*}(x)$ is the leading contribution to $t_{e}(x)$, i.e. $t_{e}^{*}(x)\approx t_{e}(x)$. Using that $t_{e}=1/\varepsilon$ we can write $\Delta E(x)=E(\varepsilon(x^{\prime}))-E(\varepsilon(x))$ as

[TABLE]

Introducing $\delta_{i}(x)\equiv\varepsilon^{\prime}-\varepsilon$ in the equation above and solving for it we find

[TABLE]

To guarantee that the two states are within $\delta_{i}(x)$ at time $t_{i}$, the distance $\delta_{x}(x)$ between $x^{\prime}$ and $x$ at time $t=0$ should be (according to Eq. (8))

[TABLE]

The result above depends on $t_{i}(x)$, while we would like to express it as a function of our observable $t_{e}(x)$. After injection the divergence is very small, and thus we replace $\lambda_{t_{i}}t_{i}$ by $\lambda_{t_{e}}t_{e}$ but we additionally multiply this expression by the total divergence $D(\varepsilon)\propto t_{e}(\varepsilon)^{2}$ discussed above (the total divergence $D(x,t_{e})$ is the product of the divergence up to $t_{i}$ and the divergence from $t_{i}$ to $t_{e}$). Our final results is then

[TABLE]

which relates the scale that two states $x$ and $x^{\prime}$ should be in order to achieve a given expected variation in the observable $E$. This expression replaces Eq. (8) for the Pomeau-Maneville map (12) considered here. The essential new feature is the appearance of the multiplying term $t_{e}(x)$. The term $e^{-\Delta E_{x}}-1$, which ideally should incorporate the desired $\Delta E$ computed using Eq. (7), does not show a strong dependence on $t_{e}(x)$ so that in practice we considered it to be a constant (i.e., we incorporate it in the proportionality constant $\Delta$).

We test the accuracy and usefulness of the results above through numerical simulations of a flat-histogram (3) simulation, as reported in Fig. 4. The success on the estimation of $P(t_{e})$ can be seen on: (i) the agreement with the traditional uniform sampling, not only in the tail of $P(t_{e})$ but also for short times; and (ii) the polynomial scaling of the computational efficiency (round-trip scales as $E^{2}$ for increasing $E$). This confirms that the proposal distribution derived in Eq. (19) achieves its goal in obtaining an efficient Metropolis-Hasting simulation. While the specific derivation presented here is valid for the Pomeau-Maneville map only, for which analytical results exist, the reasoning of deriving a $\delta_{x}(x)$ that allows to change $\log t_{e}$ can be applied more generally to maps with marginally unstable points. The main conclusion here is that the methodology summarized in Sec. II is applicable to weakly chaotic open systems with power-law distributions of the escape time, and reinforces the thesis that the strength of this methodology lies in its ability to construct the proposal distribution from the properties of the dynamical system under investigation.

III.3 Standard Map

The previous section considered an example of a weakly chaotic system for which, due to its simplicity, it was possible to derive a relationship between the distance of two trajectories $x^{\prime}-x$ and the difference in their respective observables, $E_{x^{\prime}}-E_{x}$. This motivates us to consider a more challenging problem, which shows a similar type of survival probability distribution $P(t_{e})$ but for which no analytical results are known. We consider the area-preserving standard map $(p_{t+1},\theta_{t+1})=F(p_{t},\theta_{t})$ given by OttBook

[TABLE]

This map is a paradigmatic example of the KAM scenario of mixed phase-space Hamiltonian systems OttBook ; Zaslavsky2002 . For the parameters $K$ we use in our numerical investigations ($K\gtrapprox 2$), chaotic and regular components with non-zero measure coexist in the phase space. Introducing an exit region in the chaotic component, the escape time distribution of trajectories started in the same component follow $P(t_{e})\sim t_{e}^{-\alpha}$, as shown in Fig. 5 (which again motivates the observable $E_{x}=\log t_{e}(x)$). The intermittency (or stickiness) in the standard map is weaker than the one observed in the Pomeau-Maneville map: there is (normalizable) invariant measure $\mu$ (the phase space area), the exponent is $2<\alpha<3$, and the average escape time $\langle t_{e}\rangle$ exists. The existence of a universal asymptotic exponent $\alpha$ has been long conjectured Karney ; Chirikov and the subject of extensive theoretical and numerical investigations Meiss1986 ; Zaslavsky2002 ; Cristadoro2008 . The numerical investigations in these works considered uniformly chosen initial conditions and our goal is to investigate whether more efficient choices of initial conditions can be obtained through our approach.

We start repeating the steps performed for the two previous maps to investigate how the distance between $x$ and $x^{\prime}$ changes with $t_{\star}$ in Eq. (14). The results shown in Fig. 6 show that, contrary to what was observed in the previous maps, trajectories remain close to $x_{t}$ at least up to a time $t_{\star}$, as designed in original formulation (8). Comparing to the results in Fig. 3 for the Pomeau-Maneville map, there seems to be no special time $t_{i}$ at which trajectories behave fundamentally different.

The result above indicates that $\delta_{x}(x)$ in Eq. (8) can be used to control the time up to which trajectories are close. Even without a theoretical result indicating what $t_{\star}(x)$ should be used, we see from the results above that $t_{\star}=0.5t_{e}$ is an heuristic choice that guarantees $x^{\prime}$ with $E_{x^{\prime}}$ both higher and lower than $E_{x}$ for different values of $t_{e}(x)$. Implementing this into a Metropolis-Hastings sampling method we systematically observe that the acceptance of our method approaches zero, ruining the efficiency of our sampling algorithm. To understand why this happens, recall that that a crucial simplification in the derivation of Eq. (7) was that of reversible proposal distribution g(x’—x), Eq. (4). It will be argued below that the proposal with $t_{\star}(x)=0.5t_{e}(x)$ in the standard map guarantees a bounded $\pi(E_{x^{\prime}})/\pi(E_{x})$, but it fails to guarantee a constant acceptance ratio, Eq. (1), because of the mismatch in the proposal distributions $g(x^{\prime}|x)$ and $g(x|x^{\prime})$, in violation of assumption (4). In fact, results in Fig. 6 suggest that, even though the escape time for $t_{\star}=0.5t_{e}$ lead to similar $\log t_{e}$, their respective FTLE varies dramatically. For a local proposal drawn from a half normal distribution – see Appendix A – the ratio of the proposals $g$ is given by

[TABLE]

By definition, $x^{\prime}$ is constructed to be drawn such that $|x^{\prime}-x|\approx\delta_{x}(x)$. Thus, the ratio $g(x|x^{\prime})/g(x^{\prime}|x)$ essentially depends on the ratio $R(x,x^{\prime})\equiv\delta_{x}(x)/\delta_{x}(x^{\prime})$. This allows to write $g(x|x^{\prime})/g(x^{\prime}|x)=f(R(x,x^{\prime}))$ as

[TABLE]

This function fulfills $f(1)=1$ and $f(0)=0$, and decreases to zero for $R\ll 1$ and $R\gg 1$. Thus, the more the distributions differ, the larger/smaller $f(R)$ is (depending on whether $R<1$ or $R>1$). To guarantee a constant acceptance, the proposal distribution also needs to guarantee a bounded ratio $g(x|x^{\prime})/g(x^{\prime}|x)$, which thus equates to guarantee a bounded $R(x^{\prime},x)$. Since there are no more free parameters of the proposal distribution, what remains to be analyzed is whether $R(x^{\prime},x)$ is bounded or not. Figure 7 illustrates the different values of $R(x^{\prime},x)$ obtained from the same points used to construct the histogram of Fig. 6, for the case $t_{e}=16437$. As anticipated, it indicates, that the ratio $R(x^{\prime},x)$ is orders of magnitude different from 1 specially with $t_{e}(x^{\prime})>t_{e}(x)$, which, from the preceding discussion, leads to an arbitrarily small acceptance rate.

In summary, this section showed how the proposal distribution with $t_{\star}(x)=0.5t_{e}(x)$ in the open standard map allows to propose states with an increasing $\log t_{e}(x^{\prime})$. This can be used, e.g., for algorithms that aim to find long living trajectories. However, this proposal distribution leads to vanishing acceptance rate in a Metropolis-Hastings simulation, which implies that it is not suitable to sample long living trajectories. To understand the reason for the vanishing acceptance, consider a proposed move from $x$ to $x^{\prime}$ that leads to the desired local increase of $E=\log t_{e}$. In order for this move to be accepted, as discussed in Sec. II, the probability of the reverse move (from $x^{\prime}$ to $x$) should not be vanishingly small (ideally, $g(x|x^{\prime})\approx g(x^{\prime}|x)$). However, a local move in $E=\log t_{e}$ is a large move in $t_{e}$ (for large $E$), which is multiplied by the FTLE $\lambda_{t_{\star}}(x)$ and exponentiated to compute the characteristic search scale $\delta_{x}(x)$ in Eq. (8). In the standard map, the FTLE is not sufficiently small to compensate for this increase and therefore the value of $\delta_{x^{\prime}}(x^{\prime})\ll\delta_{x}(x)$, leading to $g(x|x^{\prime})\ll g(x^{\prime}|x)$. The difference to the Pomeau-Maneville map is that the stickiness in the standard map is weaker ($\alpha>2$, larger FTLE) and more complicated Meiss1986 ; Cristadoro2008 (e.g., not a single trapping point). A possible strategy to obtain an efficient sample is to include the ratio of the proposal distributions on the acceptance rate, potentially leading to an extension of Eq. (7) which sets a tighter condition for a bounded acceptance rate.

IV Conclusions

Sampling rare trajectories in chaotic systems invoslves coupling two dynamical systems: the deterministic system we aim to study and the stochastic sampling method we construct. In a Metropolis-Hastings Monte Carlo sampling, the efficiency of the sampling depends critically on how the coupling is set through the choice of the proposal distribution. The ideas presented in this paper show how to construct such proposal distribution based on the properties of the deterministic system. In particular, we discussed how weakly-chaotic properties of the dynamics affect the proposal distribution under naive assumptions of strong chaos.

The main computational problem we discussed was to sample trajectories with very large escape time in open (weakly-chaotic) systems. For the case of one-dimensional maps with marginal points, we were able to obtain an efficient Metropolis-Hasting method to sample trajectories. For the case of area-preserving maps with mixed phase space, we showed how our approach is able to find long-living trajectories in the system but that the Metropolis-Hasting method fails due to a low-acceptance ratio. We hope these results will trigger further work on the application of Monte Carlo methods in deterministic dynamical systems, in particular to the open problem of having an efficient sampling method to estimate the tails of the survival probability in intermittent systems (e.g., Hamiltonian systems with mixed phase space in arbitrary dimensions).

V Acknowledgments

The results of this paper were obtained while JCL was a PhD student jorgeThesis at the MPIPKS in Dresden, funded by Erasmus Grant No. 29233-IC-1-2007-1-PT-ERASMUS-EUCX-1. EGA was funded by the University of Sydney bridging Grant G199768.

Appendix A Half-Gaussian local proposal

Our local proposal consists in perturbing $x$ by a finite amount $\vec{\delta}$, $x^{\prime}=x+\vec{\delta}$, characterized by a direction $\hat{\delta}$ and a norm $\delta$, $\vec{\delta}\equiv\hat{\delta}\delta$, $g(x^{\prime}|x)=x+P(\vec{\delta}|x)$. A common case is when the probability distribution is separated in two independent terms RobertCasellaBook :

[TABLE]

$P(\hat{\delta}|x)$ is uniformly distributed in the $D$ directions, and $P(\delta|x)$ has zero mean (i.e. an isotropic proposal). Additionally, we consider that $P(\delta|x)$ is characterized by a well defined scale, e.g. a half-normal distribution with mean $\delta_{x}(x)$:

[TABLE]

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