Improvements to the use of the Trajectory-Adaptive Multilevel Sampling algorithm for the study of rare events

The Trajectory-Adaptive Multilevel Sampling (TAMS) is a promising method to determine probabilities of noise induced transition in multi-stable high-dimensional dynamical systems. In this paper, we focus on two improvements of the current algorithm related to (i) the choice of the target set and (ii) the formulation of the score function. In particular, we use confidence ellipsoids determined from linearized dynamics in the choice of the target set. Furthermore, we define a score function based on empirical transition paths computed at relatively high noise levels. The suggested new TAMS method is 5 applied to two typical problems illustrating the benefits of the modifications.

especially near phase transitions, the variance of the estimated probability can peak and the convergence of the algorithm can be slow.
The aim of this work is to propose improvements to the use of the TAMS algorithm to be able to compute transitions in multi-stable systems more efficiently. The first type of improvement is the choice of the target set, which is often determined from rather arbitrary thresholds. This choice also raises more broadly the question of a precise definition of what we consider a 40 noise-induced transition between two (stable) states. In the second type of improvement, we propose a more systematic method of defining a score function, i.e., based on empirical transition paths. The modified TAMS method is first applied to an idealized gradient system, and then to a system representing a box model of the AMOC .
In section 2, we describe the methods developed to improve the TAMS algorithm. In section 3 we show how to incorporate these techniques into the definition of the score function and present the results for idealised dynamical systems and the AMOC 45 model. A discussion follows in section 4, assessing the strengths and the limitations of the new TAMS method.

Transition probabilities using TAMS
We consider finite-dimensional dynamical systems described by stochastic differential equations (SDEs), of the following form: where X t ∈ R n and F : R n → R n is the drift field. The noise term W t ∈ R m consists of m independent Wiener processes with the matrix G : R n × R m being the noise matrix.
A prominent example of such a system is a model of a free particle moving in a two-dimensional double-well potential (with n = m = 2). The drift term in the time-evolution equation for the variables x and y is in this case where V (x, y) represents the potential ( Fig. 1(a)). In the deterministic case (i.e. G = 0), the stable steady states of the system are X A = (−1, 0) and X B = (+1, 0), while the unstable steady state is X C = (0, 0). Without fluctuations, a particle starting in X A stays in X A . However, because of the presence of unresolved processes (such as thermal fluctuations), modelled by the noise term in the SDE, the particle can move away from X A and in some cases make the transition to the state X B . An example 60 of such transition is shown in Fig. 1(b). The transition probability that a trajectory starting in X A reaches a neighbourhood B around X B before time T is indicated by P (τ B < T |X 0 = X A ). Here, τ B denotes the stopping time associated with reaching B. The TAMS algorithm is based on a selection and mutation process of discarding and branching trajectories, which are ranked according to a score function φ : R n × R → R. For a given state (X, t) ∈ R n × R, φ(X, t) is supposed to measure how likely it is to start in (X, t) and reach B before time T . As a consequence, if the choice of score function is successful, the probability that a trajectory reaches the target set B keeps increasing at each iteration of the algorithm, which is why this method is more efficient than brute-force techniques. A visual representation of the algorithm is given in Fig. 2

and a step-by-step description is provided in Appendix
A (more details can be found in Lestang et al. (2018)). trajectories are ranked according to their score φi, which is the maximum value of the score function φ along the trajectory. Then, at each iteration, the trajectory with lowest score (1 in the figure, with score φ1) is discarded. It is replaced by picking a trajectory uniformly at random from the other ones (2 in the figure), then computing the earliest position at which it reaches a higher score than the discarded trajectory (orange dot) and finally using this position as the branching point for the new trajectory (1 ). This is repeated until all the trajectories reach B or the number of iterations reaches a predefined limit.

70
Consider a general SDE (1) and two stable states X A and X B of the corresponding deterministic system. The TAMS algorithm (cf. Appendix A) needs a score function to be defined, which allows to rank trajectories, and select the ones to discard. The optimal score function φ com (X, t), i.e. the score function that minimises the variance of the probability estimator, is called the static committor. Its generic expression is given by the following conditional probability (Lestang et al., 2018):

75
where T is again the fixed duration of the trajectories and τ B is the stopping time associated with reaching the target set B. In other words, φ com (X, t) is the probability that a trajectory starting in X at time t reaches B before time T . This expression is quite natural because ideally, the score of (Y, s) should be higher than the score of (X, t), i.e., φ(X, t) ≤ φ(Y, s), if and only if . This condition is clearly satisfied by φ com (and in fact any increasing function of φ com ).
However, the expression given in eq. (3) is generally unusable because it is the very quantity that we want to compute. For 80 instance, φ com (X A , 0) is precisely the transition probability that TAMS estimates. As a conditional probability of the form P (Y, s|X, t), the committor φ com (X, t) satisfies the backward Fokker Planck equation (Lestang et al., 2018): However, solving the backward Fokker Plank equation in systems with many degrees of freedom is computationally infeasible.
Moreover, even if the committor is available on a grid used for the discretization of (4), using interpolation to evaluate it during 85 a TAMS loop can also have a prohibitive computational cost. Bouchet et al. (2019) made use of a score function based on the distances of the state X from the starting state X A and the destination equilibrium X B , respectively: it is defined as Alternatively, a Gaussian-shaped score function φ gauss was proposed in Baars (2019). It is defined as: Here β ∈ R is a parameter controlling the decay and X C is the saddle state of the system.
Due to the general expressions of φ gauss and φ dist , these score functions can be used in systems of any dimension.

Definition of the target set
Once the score function has been chosen, a threshold needs to be defined for the TAMS algorithm to converge, so that the 95 occurrence of a transition can be detected. In other words, we do not expect each trajectory that undergoes a transition to reach exactly the destination equilibrium X B , but rather a neighbourhood of it (B). The target set B can then be defined according to a level set φ target of the score function φ: However, different score functions φ and different level sets φ target correspond to different target sets B, which can differ in 100 volume and in shape. Moreover, often the level set φ target is defined somewhat arbitrarily. For example, by using φ gauss with target score φ target = 0.85, 0.9 or 0.95, we found that the average of the transition probability estimator for a two-dimensional double-well potential system (2) can vary up to 30%. This may not be a concern if one only cares about the order of magnitude of the transition probability but it will be problematic if quantitative comparisons are needed. Moreover, a poor choice of target set can lead to inaccuracies when a trajectory has a score greater than φ target without actually making the transition in all 105 the degrees of freedom. Thus, there is a need for defining a canonical choice of a target set B, which needs to fundamentally address what is considered to be a noise-induced transition.
For this purpose, we use the concept of confidence ellipsoid. This is an ellipsoidal neighbourhood around a stable equilibrium state, inside which a trajectory subject to the locally linearized dynamics stays, within a certain confidence level (Cowan, 1998). Consider a general SDE system given by (1). Because the drift vanishes at an equilibrium state F (X B ) = 0, its first 110 order approximation around the equilibrium state X B via Taylor expansion is: where A(X B ) = ∇F (X B ) ∈ R n×n is the Jacobian matrix of F at X B . Using a translation X = X − X B , the first order approximation of the SDE is then:

115
Because the drift term has been linearized, this is the equation for an n-dimensional Ornstein-Uhlenbeck process. The stationary probability density function (PDF) f of the approximating process is Gaussian (Cowan, 1998) and given by where C B ∈ R n×n is the covariance matrix of the system calculated in X B and . C −1 B the norm induced by its inverse The covariance matrix C B can be thought heuristically as the matrix containing the 120 correlations E(x i x j ) (with X = (x 1 , . . . , x n )), which generalises the notion of variance in n-dimensions. C B is obtained by solving the Lyapunov equation (see Kuehn, 2012, for the full derivation): We then define the confidence ellipsoid, which has C −1 B as shape matrix, as follows where Q α is the quantile of confidence level 1 − α of the n-dimensional χ 2 distribution (Cowan, 1998); usually 1 − α = 0.95.
The directions of symmetry of the ellipsoids are given by the eigenvectors of the covariance matrix C B and the radii are given by the corresponding eigenvalues and the confidence level 1 − α. Intuitively, the greater the eigenvalue, the more a trajectory fluctuates in the given direction.
The (1 − α)-covariance ellipsoid represents the n-dimensional volume where a trajectory is confined with confidence level 1 − α, provided its dynamics is well approximated by the first order expansion at the equilibrium point. An illustration of the 0.95-confidence ellipsoid for the double-well potential is shown in Fig. 3. The confidence ellipsoid E constitutes a way to meaningfully define the target set B with minimal arbitrary parameters. As shown in Fig. 3, when initializing a trajectory at X A or X B in the two-dimensional double-well system, it stays inside the correspondent ellipsoid with a certain confidence 1 − α = 0.95. In the next section, we show how to incorporate this choice of target set in the score function φ.

Estimating the typical transition path using histograms
The second line of improvement of the score function concerns the estimation of typical transition paths of the dynamical system. In the zero noise limit, the Freidlin-Wentzell theory of large deviations predicts that transition paths cluster around the most probable transition path, called the instanton (Freidlin and Wentzell, 1984). On the other hand, in the finite noise regime, transition paths may deviate from the instanton. Moreover, instantons may be computationally inaccessible for systems with 140 many degrees of freedom. Therefore, it can prove more relevant to estimate empirically the typical transition path that the system follows at a given finite noise level, which is the approach we follow here.
The idea is to first accumulate transition paths at a noise level where transitions are frequent enough (typically p > 10 −3 ) so that any sampling method (direct Monte Carlo or TAMS with naive score functions) can be used. Then, we compute the spatial histogram of the transition paths over a discretized phase space using n−dimensional boxes. This provides the 145 spatial distribution of the transition paths, which is concentrated around a typical transition path, reminiscent of an instanton phenomenology, which was also observed in more complex systems (Bouchet et al., 2019). From the spatial histogram, we extract a typical transition path. The main steps of the path-finding algorithm are listed below: (i) the trajectory of the typical transition path starts in the box of the histogram containing the initial state X A ; (ii) the next box in this trajectory corresponds to the neighbour which has the highest nonzero histogram value but which has 150 not already been visited by the typical transition path; (iii) the algorithm stops if it reaches the box containing the target state X B .
In addition, the full typical path estimation algorithm uses a self-correcting method to avoid dead ends when there are no valid neighbours to be the next point in the trajectory. The spirit of the path-finding algorithm is similar to the depth-first search algorithm (Cormen et al., 2009). We found that, as long as the histogram is not fragmented, i.e., there is a sufficient number of 155 accumulated trajectories or large enough histogram bins, the algorithm converges. Possible artefacts created by this estimation include spiralling near the initial equilibrium (because of the concentric shell structure of the local probability density function (10) and zigzagging at the histogram box size. They can be both addressed by a clean-up algorithm: starting from the first box, at each box X j , if the trajectory goes back to one of its neighbours at a later time, with X j+k being the latest neighbour visit, we erase the points X j+1 , . . . , X j+k from the trajectory.

Results
In this section, we apply both modifications to TAMS (ellipsoids in the score function and typical path estimation) to different problems.

Incorporating ellipsoids in the score function
First, we apply the modified TAMS method to the two-dimensional double-well system defined by (1). Here we see, which will 165 be more general, that the level sets of the score function do not have the shape of an ellipsoid. Hence, there is no level φ target Here, we propose a general method to modify any target score function so that we are able to choose the target set to be exactly the confidence ellipsoid of X B . Let to the procedure described in eq. 14. The level set φ target (dotted blue line) of the former score function φgauss is tangent to the ellipsoid.
The level set φ target (dotted red) of the modified score function φgauss coincides with the ellipsoid. (b) Same plot, for the score function E be the 0.95-confidence ellipsoid around the equilibrium state X B and φ (e.g. φ gauss ) the score function to be modified. We first compute the level φ target , defined as the minimum of the score function φ on the confidence ellipsoid E: such that the set {X ∈ R n | φ(X) > φ target } contains the ellipsoid E and is tangent to E. This can be done numerically by generating a mesh of points around X B , then selecting the points inside E by comparing their norm with the quantile Q α and finally computing the minimum φ target of φ on these points. Then, define the modified score function φ in the following way: The target set B = {X ∈ R n | φ(X) > φ target } for the modified score function φ turns out to coincide with E. We apply this procedure on both the score functions φ gauss and φ dist and the results for the improved score functions φ dist and φ gauss are shown in Fig. 4.

Designing a score function based on a typical transition path 180
In order to show how to design a score function based on a typical transition path, we consider a two-dimensional system slightly less trivial than the double-well system, i.e. a two-dimensional system with the following potential: depicted in Fig. 5(a). It consists in two energy minima at X A ≈ (−5.77, 0) and X B ≈ (5.77, 0) and a potential barrier spanning y < 15 and at x = 0. The dynamics is then given by the SDE (1) The dynamics of the system is quite interesting, as it exhibits two distinct regimes for transition paths, depending on the noise level σ (assuming G = σI). At high noise (σ 2 ∆V , where the potential barrier height is  Rolland and Simonnet (2015) investigated the convergence properties of AMS using a triple-well potential system, 190 which also exhibits two regimes of preferred transition paths depending on the noise. They found a strong dependency of the statistics of the rare-event algorithm (e.g., the number of iterations) and the duration of reactive trajectories on the choice of score function. We expect the same behaviour when applying TAMS to the system with potential (15) and that this system reveals differences in performance between various score functions. already mentioned, some artefacts created by the estimation (such as spiraling or zigzagging) can be corrected using a cleanup algorithm. The result is shown in Fig. 6(b): the empirical estimation of the typical path resembles the instanton around which the transition paths are clustered at lower noise ( Fig. 6(b)). The instanton was calculated by implementing the geometric minimum action method (Heymann and Vanden-Eijnden, 2008). Note that if, instead, we accumulated trajectories at high noise 200 σ > 10, we would obtain trajectories going from X A to X B in a straight line, which is the typical path at high noise similar to Rolland and Simonnet (2015). This typical transition path is then substantially different from the instanton, which goes through the upper channel. Thus, our method can be advantageous in multistable systems where the typical transition path depends on the noise level. We can start at high noise and reapply the empirical estimation of the typical path each time the noise level is decreased.

205
Given a typical transition path C, we present the design of a score function φ C which encourages trajectories to follow the transition path C such that it gives a reasonable approximation of the static committor. Let us consider a trajectory C(s) in R n parametrised by arclength s ∈ [0, 1]. Then we define the score function φ C , called the path-based score function, such that it grows from 0 to 1 along the trajectory from X A to X B and decays exponentially along the direction transverse to the trajectory: (free parameter). The score function φ C is shown in Fig. 7 for the estimated transition path C shown in Fig. 6 and two values of decay length d 0 = 20 and d 0 = 200. The score function increases from 0 to 1 along the trajectory. Thus, it encodes the 215 preferred direction that the system has to follow. This contrasts with the generic score functions φ gauss and φ dist which are symmetrical in y: they do not contain the information that the system has to increase in y in order to make the transition to X B . Note that the method we developed here can be applied, in principle, to systems of any dimension. As shown in Fig. 7(b) the score function φ C is discontinuous because the trajectory has positive curvature. The discontinuity is located near the axis x = 0. Indeed, when crossing the axis x = 0, the closest point on C changes and s(X, C) is discontinuous.

220
Next, we applied TAMS with the path-based score function φ C to the two-dimensional system (15). We compare its performance with the previously defined score functions φ dist and φ gauss . In fact, we use the associated modified score functions, such that the target set B matches the 0.95-confidence ellipsoid (we drop the tildes for readability). We use the following parameters: duration of a trajectory T = 20, time step dt = 0.01. We show in Fig. 8 The performances of the numerical methods are next measured using the work-normalised relative error , which combines the variance of the algorithm and its computational cost (Glynn et al., 2009):

230
where µp and σp are the mean and standard deviation of the probability estimate over the different instances, ω is the average number of time steps calculated in one realisation. In short, measures how precise the numerical method is at equal computational cost. The smaller , the better the algorithm performs.
The results are shown in Fig. 8(b). For the system (15), the score functions φ dist , φ gauss and φ C have little difference in performance. For the lowest noise values, the path-based score function has at most a 30% smaller error than the score 235 function φ gauss . Changing the decay length d 0 hardly changes the error . When adjusting the parameters of the potential V or applying the same method to the triple well system used in Rolland and Simonnet (2015), the performance gain, while being often present, never systematically exceeded 30%. All in all, in this category of two-dimensional systems, using the path-based score function approach yields little improvement. Transition probability p Work-norm. rel. error (as shown in Appendix B). This model can be formulated as

245
In the equations above, we split the state of the system X t into two parts: Y t , which includes salinities of four of the boxes plus the depth of the pycnocline D, and Z t , which represents the salinity of the deep box (S d ). As the noise is applied only on the asymmetric component of the atmospheric freshwater flux (E a ), it directly affects only two of the variables S n and S s ). Moreover, the stochastic increments associated with the two variables are identical, to make sure that each decrease of freshwater forcing in the southern box results in the same increase of it in the northern box in the model. Therefore, the 250 noise is not spatially independent and the noise matrix G is no longer diagonal: it consists of a (5 × 1) row vector, with only two elements different from zero. For a reasonable choice of the parameters, the deterministic system is in a bistable regime , which means that there are two possible equilibrium states under the same forcing conditions. In general, we are interested in studying transitions between the present-day AMOC (X A ) and the collapsed state (X B ).
For a differential-algebraic system of equations (DAEs), such as the system in eq. (18), we need to be particularly careful while computing the covariance ellipsoids. First of all, we make use of the Schur complement of the Jacobian of the system, which allows to calculate the covariance matrix when an algebraic constraint is present (Baars et al., 2017). Nevertheless, the resulting matrix C B is singular, with two eigenvalues being equal to zero. One of the corresponding eigenvectors is a vector pointing in the direction of the variable D (depth of the pycnocline). The reason behind it is that the differential equation governing the evolution of D does not contain any of the other variables (see Castellana et al., 2019), when the system is in the 260 collapsed state (X B ). This results in D not being affected by the noise, as this is imposed only on two of the salinities. Hence, we compute the covariance matrix C S B relative to the salinities, removing one degree of freedom from the original matrix. Unfortunately, such matrix still has one zero eigenvalue, which is due to the salinity conservation (the algebraic equation in the system (18)). To overcome this problem, we compute the Moore-Penrose inverse (or pseudoinverse) of the covariance matrix,  Fig. 9.
For the system (18), the modified score function is more complicated, as the covariance matrix used to construct the ellipsoid contains only the degrees of freedom related to the salinities of the model, leaving the variable D (depth of the pycnocline) out. From a geometric point of view, that means that the covariance ellipsoid around X B is degenerated along the D-direction. 270 Fig. 9 shows a projection of the ellipsoid -once the noise amplitude is fixed -on the plane identified by the variables S n and S s (respectively, the salinity of the northern box and the one of the southern box in the model). The projection was obtained calculating the conditional covariance matrix of the two variables into consideration, given that the other variables are set on their mean value (Wasserman, 2013). The two-dimensional ellipse contains the large majority of the projected time points (on the same plane) of a trajectory that wanders around the equilibrium. By construction, the confidence level of the confinement is 275 higher than the one prescribed for the full dimensional ellipsoid (in this case 0.95). Clearly the level sets of the score function φ gauss do not coincide with the one of the ellipsoid. Hence, the importance of modifying the score function appears evident.
When constructing an improved score function, in order to evaluate if a state belongs to the neighbourhood of X B , we need to check two conditions: (i) whether the state of the system is inside the salinity covariance ellipsoid drawn around the destination equilibrium, and (ii) we need to verify that the variable D of the state is the same as the one of X B . Hence, the improved score 280 function for the box model reads where X S indicates the part of the state vector representing the set of the salinities, whereas X D represents the variable D. As already mentioned, to check whether a certain state belongs to the salinity ellipsoid, we made use of the pseudo-inverse of C S B in the definition (12). To be able to assess the relevance of a proper definition of the target set in the TAMS algorithm, and hence the importance of using the improved version of the score function, we computed transition probabilities of the AMOC from the present-climate state to the collapsed state for reasonable values of the atmospheric forcing and noise. In particular, we choseĒ a = 0.20 Sv and f σ = 0.16. This last value represents the ratio between the standard deviation of the noise in the atmospheric forcing and its mean value . The number of trajectories used in the algorithm was set to 100, and the time scale at 290 which the probabilities were evaluated was chosen to be 1, 000 and 10, 000 years, respectively. For each transition probability, we used three different versions of the algorithm, based on three different settings for the score function. The first two versions were implemented using φ gauss in (6), with two different choices for the threshold φ target . In the third version, we used φ box as defined in (19), where the starting score function was φ gauss . Each probability was calculated running 15 instances of the algorithm, and then computing the median and the interquartile range (IQR). The results are shown in the Table 1 below.

295
When running TAMS to compute transition probabilities between two states of the Atlantic Circulation in this model, with different versions of the score function, we found a considerable discrepancy between the obtained values. In particular, it appears that setting a very high threshold in the Gaussian score function makes the algorithm detect no transitions (we set up the algorithm so that it stops when the probabilities involved are smaller than 10 −9 ): the reason behind this is that, because of 300 the presence of the noise, we don't expect the state of the system to stay indefinitely close to the destination equilibrium, but rather to wander around it. Therefore, the score function, which assigns the maximum score only to a very small neighbourhood of the equilibrium, is not able to properly recognise transitions. Moreover, it is not surprising that, when using φ gauss with a Score function p at 1000 years [IQR] p at 10000 years [IQR] φgauss, φ target = 1 − 10 −4 < 10 −9 < 10 −9 φgauss, φ target = 1 − 10 −2 (1.1 [0.6 : 1.2]) × 10 −3 (5.3 [4.7 : 6.0]) × 10 −2 φbox < 10 −9 (2.0 [0.7 : 3.0]) × 10 −3 Table 1. Results of the transition probabilities for the AMOC model using different score functions smaller value of φ target (0.99) or φ box , we obtain different results: as the shape of the covariance ellipsoid is not spherical (see Fig. 9), we expect φ gauss to detect transitions even though the state is actually still far from the destination equilibrium, at 305 least in certain directions. As a general rule, we expect φ gauss to give results different from φ box as long as the ellipsoid of the system is not spherical, regardless of the choice of the threshold φ target .

Summary and discussion
We presented and applied several improvements to the TAMS rare-event algorithm, when used to compute transitions in multistable systems. The first improvement was based on a more rigorous criterion to define noise-induced transitions involving 310 confidence ellipsoids E to formalise this criterion. In turn, this led to the rigorous choice of the target set B = E B which was traditionally set by rather arbitrary thresholds. We then showed how to incorporate this definition of B into the score function φ. For certain classes of systems, like the ones containing algebraic constraints in addition to differential equations, or when the noise does not affect one or more directions in the variable space (i.e. the associated covariance matrix is singular), this method requires some precautions. In particular, for the box model in Castellana et al. (2019), we had to adapt the definition of 315 the improved score function (19), as well as calculate the pseudo-inverse matrix of the covariance matrix, in order to compute the ellipsoid.
This method, while being quite general, has several limitations. While the modified score function φ is continuous outside Fig. 4). This means that in the TAMS algorithm, branching will never occur inside M , but at the boundary ∂M . This can have an influence on the convergence of TAMS if the 320 level sets of the initial score function φ have a pathological shape near X B and the spatial extension of D is not negligible.
Nevertheless, we expect this to have little impact because this phenomenon is localised near the target state X B . Therefore, the trajectories will naturally converge towards X B as a result of the dynamics, even without the help of the branching process of TAMS. However, to ensure that the confidence ellipsoid E defines a meaningful target set, one needs to be sure that E is contained inside the basin of attraction of the target state X B . While this is the case in the limit of small noise σ → 0, it might 325 not be the case for finite noise. A solution to this issue would be to compute the basin of attraction V of X B and define the target set B as the intersection B = E ∩ V. However, we expect that in the generic case, this occurs when the noise level σ is high enough so that transitions are less rare and a direct Monte Carlo estimation is sufficient to estimate transition probabilities.
Next, we proposed a systematic method of defining a score function, designed to approximate the static committor, based on empirical transition paths. We proposed an algorithm to estimate the typical transition path under a high noise level, which is 330 then used to define a family of score functions with a single decay parameter d 0 . We applied our method to a two-dimensional well with a potential barrier. We found that our typical path estimation gave satisfactory results and that the associated score function φ C , while discontinuous, remained unbiased and relatively insensitive to the change of decay parameter d 0 . While we did not find significant performance improvements over existing non-trivial score functions, we think that differences will become apparent if applied to higher dimensional systems, where there are more directions to fluctuate in. We did not show 335 any results of the application of this method to the box model of the AMOC , because for this rather complicated case, we were not able to efficiently compute typical transition paths from the trajectory histograms of the system.
One key limitation of our approach of constructing the path-based score function φ C is the computer memory needed to store the trajectory histogram, which becomes prohibitively huge for high-dimensional systems such as discretised partial differential equation (PDE) systems. As an example, a 50 × 50 two-dimensional grid storing 4 variables in each cell (e.g. two 340 velocity components, pressure and a tracer) with 10 bins of resolution in each degree of freedom would require more than 10 Petabytes of memory, which is unfeasonable. However, this limitation can be easily bypassed by defining the objects needed to run the TAMS algorithm, namely the score function φ and the target set B, in a reduced space of much fewer dimensions.
For instance, Bouchet et al. (2019) studied the dynamics of the barotropic beta-plane quasi-geostrophic equations describing Jupiter's turbulent atmosphere. While the PDE system was evolved in the full phase space, their rare-event algorithm was run 345 in a reduced 3-dimensional phase space defined by three Fourier coefficients. The target set B and the score function were defined on this reduced space. Moreover, they accumulated transition trajectories in a 3-dimensional histogram and showed their concentration around instantons. By applying the path-finding algorithm, an empirical estimation of the instanton can be made. This offers a viable alternative to solving a minimization problem in the full space to compute the instanton and then project it in the reduced space, which is next to intractable for this system.

350
Another way to define the reduced space V in which to run TAMS is to consider the principal components, also called empirical orthogonal functions (EOF), which are the eigenvectors of the covariance matrix. One idea, suggested by Baars (2019), is to retain the EOFs {Y A 1 , . . . , Y A k , Y B 1 , . . . , Y B k } with largest variance (i.e. eigenvalue of the covariance matrix) at the initial state X A and target state X B . EOFs represent the directions in which the system fluctuates the most. They are then assumed to be the directions which capture best the noise-driven dynamics. When studying transitions in a PDE model of the Atlantic 355 Meridional Overturning Circulation (den Toom et al., 2011;Baars et al., 2019) projected the dynamics in a reduced space W ≡ Span{X A , X B , Y A 1 , . . . , Y A k , Y B 1 , . . . , Y B k } of dimension (∼ 500) still too large to apply the histogram method directly. However, one idea is that the TAMS algorithm could be run in an even smaller space V ≡ Span{X A , X B , Y A 1 , . . . , Y A d , Y B 1 , . . . , Y B d } (of dimension <10), while still computing the dynamics in the space W. Then, the memory required to store a histogram becomes reasonable and our histogram method can be applied.

360
Another potential issue of our modified TAMS method is the fact that the score function φ C is discontinuous because the trajectory has positive curvature, as shown in Fig. 7(b). The discontinuity is located near the axis x = 0. Indeed, when crossing the axis x = 0, the closest point on C changes and s(X, C) is discontinuous. In fact, in the mathematical proofs about the where the function θ(x) is the Heaviside step function. The transports depend on the variables via the following relations: where the density of the generic box i is defined as

385
The volumes depend, in turn, on the variable D: The reference parameter values are shown in table B1.