Sources and Sinks of Rare Trajectories in 2-Dimensional Velocity Fields Identified by Importance Sampling
SSOURCES AND SINKS OF RARE TRAJECTORIES IN2-DIMENSIONAL VELOCITY FIELDS IDENTIFIED BY IMPORTANCESAMPLING
MEAGAN CARNEY*, HOLGER KANTZ † Abstract.
We use importance sampling in a redefined way to highlight and investigaterare events in the form of trajectories trapped inside a target almost invariant set. Wetake a transfer operator approach to finding almost invariant sets of a reconstructed 2-dimensional flow of the atmosphere from wind velocity fields provided by the PortableUniversity Model of the Atmosphere. Motivated by extreme value theory, we consideran observable φ ( x ) = − log( d ( x, γ )) maximized at the center γ of a chosen target almostinvariant set. We illustrate that importance sampling using this observable provides anenriched data set of trajectories that experience a rare event, in this case defined as endingnear γ . Backwards reconstruction of these trajectories provides valuable information oninitial conditions and most likely paths a trajectory will take toward a rare event.In this specific setting, one can think of an almost invariant set as a region in theatmosphere where only a small number of particles move in and out of the region. Forshort time intervals, these regions are represented physically by individual eddies in theatmosphere; while longer time intervals suggest formation by eddy interaction such asparticles trapped between the spin of two eddies. Particle movement in and interactionwith these sets can provide useful information on the most probable path of a storm (inthe former case), or origin of pollution (in the latter), for example. Introduction
Atmospheric eddies play a major role in extreme weather phenomenon such as heat-waves, hurricane movement, and pollution distribution [15, 16]. Topologically, these eddiescan be seen as almost invariant sets of a flow where there is minimal movement of particlesfrom within and outside of a region [8]. Understanding where these eddies occur and thelikelihood of trajectories ending up inside them can provide a new and useful perspectiveon atmospheric movement.We reconstruct a 2-dimensional model of atmospheric flow defined on a space Y ⊆ R from wind velocity fields provided by the Portable University Model of the Atmosphere(PUMA) [13]. Following recent literature, we estimate the transition probability matrix Date : June 16, 2020.*Max Planck Institute for the Physics of Complex Systems, N¨othnitzer Str. 38, D 01187 Dresden,Germany. Email: [email protected] † Max Planck Institute for the Physics of Complex Systems, N¨othnitzer Str. 38, D 01187 Dresden,Germany. Email: [email protected] thanks to Matthew Nicol for his expertise and advice on the foundations of this paper. Thanks toFrank Lunkeit for helpful discussions and information on the Portable University Model of the Atmosphere. a r X i v : . [ n li n . C D ] J un MEAGAN CARNEY*, HOLGER KANTZ † (TPM) of the flow on a subspace X ⊂ Y by taking a fine grid of boxes and measuringtransitions of particles from one box to another over a finite time interval. Intuitively, onecan think of an almost invariant set as a region where only a small number of particlesmove in and out of the region. We take a transfer operator approach to finding thesealmost invariant sets by approximating the Perron-Frobenius operator with the associatedTPM. Spectral properties of this operator provide information on the invariant structureof the space. This approach is discussed in detail in [5, 6, 4, 7] and applied numericallyto an ocean flow model in [8]. For adaptation purposes, our methods differ slightly fromthose of the listed literature by including variations on the TPM [14] and employing aspectral clustering approach equipped with K -means [17, 1, 2]; however the foundationalarguments remain the same. We discuss this foundation and set of differences in detail inthe following sections.Motivated by extreme value literature [10, 12], we consider an observable φ ( x ) = − log( d ( x, γ ))where φ : Z → R is defined for every x ∈ Z , X ⊂ Z ⊆ Y and γ is the euclideancenter of the almost invariant set. In this way trajectories of the observable under theflow are maximized as they approach the center of the almost invariant set. Under someflow f t : Y → Y , it is often of interest to consider a set of random variables defined by X N = φ ◦ f t ( x N ) for a set of initial values x N ∈ Z at some fixed time t . For our choice of f t and φ , we show that ( X N ) behaves as though it comes from some shifted exponentialdistribution where f t ( x N ) → γ gives X N → ∞ .We apply an importance sampling method, called genealogical particle analysis (GPA),that exponentially tilts the distribution of ( X N ) so that the probability of observing largervalues (and hence, values of f t ( x N ) closer to γ ) is increased [3, 18, 15]. GPA works bykilling and cloning trajectories under the flow at specified sampling times based on aweight function that determines the performance of a trajectory. Large values of theweight function indicate that a trajectory behaves as though it comes from the target(tilted) distribution. In the end, the surviving trajectories represent the set that has ahigher probability of ending near γ . Backwards reconstruction of these trajectories allowsus to find the set of most likely paths that end in the almost invariant set within a specifiedtime interval.2. Transfer Operators and Almost Invariant Measures and Sets
We begin with a set of definitions from [6] that will be used to establish reasoning forthe methods used in this investigation. A function p : X × B → [0 ,
1] in space X withassociated σ -algebra B is said to be a stochastic transition function if1. p ( x, · ) is a probability measure ∀ x ∈ X p ( · , A ) is Lebesgue measurable ∀ A ∈ B If µ ∈ M , with M the space of all probability measures on B , satisfies µ ( A ) = (cid:90) p ( x, A ) dµ ( x ) ∀ A ∈ B OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 3 then µ is said to be an invariant measure of p . Under the assumption that the probabilitymeasure p ( x, · ) is absolutely continuous w.r.t Lebesgue m ∀ x ∈ X we have, p ( x, A ) = (cid:90) A k ( x, y ) dm ( y ) ∀ A ∈ B with transition density function k : X × X → R , k ( x, · ) ∈ L ( X, m ) and k ( x, y ) ≥ P : M → M is defined by,
P µ ( A ) = (cid:90) p ( x, A ) dµ ( x ) . If p is absolutely continuous with density function k , then P is defined on L by, P g ( y ) = (cid:90) k ( x, y ) g ( x ) dm ( x ) ∀ g ∈ L . A measure µ ∈ M is invariant if and only if it is a fixed point of P e.g. the invariantmeasures correspond to eigenfunctions of P for the eigenvalue λ = 1 ( P µ = λµ = µ ).If we let p ( x, · ) = δ f ( x ) with f : X → X a diffeomorphism on X then one can show bydefinition, µ ( A ) = (cid:90) p ( x, A ) dµ ( x ) = µ ( f − ( A ))Hence, in knowing the spectral information of the operator P for deterministic systems wecan obtain invariant measures for the diffeomorphism f . Of course, it is important in ourcase to consider a stochastic system rather than a deterministic one. In this sense, we maystart by introducing a small perturbation into the original system. The density function k associated with this ε > k ε ( x, y ) = 1 ε n m ( B ) 1 B (cid:18) ε ( y − x ) (cid:19) , x, y ∈ X where B = B (0) is an open ball of radius 1 in R n . Then, p ε ( x, A ) = (cid:90) A k ε ( f ( x ) , y ) dm ( y )and the Markov process defined by any initial probability measure µ and transition function p ε is a small random perturbation of the deterministic system f . However, many invariantmaps form more complicated structures in the phase space with regions of periodic andchaotic movement occupying closely neighboring spaces. In this case, perturbations of themap may cause the space to no longer be decomposed into invariant sets [11]. Fortunatelyfor certain small enough perturbations, these can be recovered as almost invariant sets.Informally, an almost invariant set is one where the probability of mapping from set A back to set A is [4], ρ ( A ) = (cid:82) A p ε ( x, A ) dµ ( x ) µ ( A ) ≈ . MEAGAN CARNEY*, HOLGER KANTZ † Furthermore, p ε ( x, · ) → δ f ( x ) as ε → (cid:90) A p ( x, A ) d µ ( x ) = (cid:90) A δ f ( x ) ( A ) dµ ( x ) = µ ( f − ( A ) ∩ A ) . If f ε is a differentiable, invertible time t map of a smooth flow, then P ε simplifies to P ε g ( y ) = g ( f − ε y ) | det Df ε ( f − ε y ) | where || P ε || = 1 and P ε g ≥ g ≥ P ε is often approximated bya finite dimensional Galerkin approximation based on a fine partition { B , . . . , B m } of X [4, 6, 8].Following the setup of [8] we let ( x, t, τ ) → f ε ( x, t ; τ ) represent the non-autonomous flowwhere [ t, t + τ ] is the full time interval of interest and τ is taken on the order of a day. Inthis way, f ε ( x, t ; τ ) is the end position in X of a trajectory beginning at x ∈ X , time t , andflowing for τ time. A trajectory x ( t ) := f ε ( x , t ; t ) is a solution to the non-autonomousODE dxdt = V ( x ( t ) , t ) with initial condition x ( t ) = f ε ( x , t ; 0) and vector field V . Forour purposes, we will only consider the case when X ⊆ R , however many of these resultsextend in an obvious way to flows on R n . The transition probability matrix that is formedunder the flow from t to t + τ given by, P t,τ,i,j = m ( B i ∩ f ε ( B j , t + τ ; − τ )) m ( B i )represents the approximation to P ε . Remark.
It is sometimes of interest to consider the flow over a restricted region X in thelarger space. In this setting, a point x may flow outside of the region X in t + τ timecausing the rows of the TPM to no longer sum to 1. Normalization of the rows has beensuccessfully applied in the literature [8]; however, a more natural setting for this study isallowing x to flow outside of X . To accomplish this, we introduce an ”external” tile [14]to the system represented by an additional row/column pair in P t,τ,i,j with the probabilityof transitioning into this tile equal to 1 − (cid:80) mj =1 P t,τ,i,j and transitioning out to any tilewith equal probability. Intuitively, the addition of this tile can be thought of as a point x flowing under the map until it reaches the boundary of X and then being sufficientlymixed before being injected back into X .The following nontrivial result from [4, Prop. 5.7] relates the measure ρ ( A ) to anyeigenvalue λ of the corresponding Perron-Frobenius operator P ε (and hence, P t,τ,i,j ), λµ ( A ) = ( ρ ( A ) + ρ ( X \ A ) − µ ( A ) . (1)When λ ≈ µ is close to the invariant measure ofthe system. In a similar way, if we consider the right hand side where the sets A and X \ A form a partition of the space then finding an almost invariant measure can be viewed as amaximization problem of both ρ ( A ) and ρ ( X \ A ). This approach is discussed in detail in[6]. OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 5 Graph Partitioning and Invariant Sets
Following the relationship described in (1), we may consider q partitions of the spacewhere we seek to maximize the value, ρ ( A , . . . , A q ) = 1 q q (cid:88) k =1 ρ ( A k )by varying the partitions A , . . . , A q such that A k ∩ A (cid:96) = ∅ for k (cid:54) = (cid:96) and ∪ qk =1 A k = X .For our ε -perturbed flow f ε , we may equivalently write the measure ρ of any set A as, ρ ( A ) = µ ( A ∩ f ε ( A, t + τ ; − τ )) µ ( A ) . Given a fine partition of the space X = { B , . . . , B m } and a corresponding transitionmatrix P t,τ,i,j let, p i = m ( B i ) m ( X ) , (2)where in practice m ( B i ) = area of B i if f ε is the defined flow on R , for example. In thissetting it is natural to define the probability measure [6] by, µ m ( A ) = m (cid:88) i =1 m ( A ∩ B i ) m ( B i ) p i . (3)In fact, µ m → µ strongly for small ε perturbed deterministic dynamical systems [4]. From[6, Prop. 6.4] we have that for a subset A = ∪ i ∈ I B i made up of discrete boxes B i ⊂ X which form the fine partition of the space of the flow f (cid:15) and I ⊂ { , . . . , m } of indices, ρ ( A ) ≈ (cid:80) i,j ∈ I p i P t,τ,i,j (cid:80) i ∈ I p i (4)Generally, P t,τ,i,j is not reversible because, p i P t,τ,i,j (cid:54) = p j P t,τ,j,i or equivalently, µ ( B i ∩ f ε ( B j , t + τ ; − τ )) (cid:54) = µ ( B j ∩ f ε ( B i , t + τ, − τ )) . In other words, the probability of being in box B i at time t and then B j at time t + τ is notequal to the probability of being in box B j at time t and then B i at time t + τ . However,if we define the time-reversed quantity,ˆ ρ ( A ) = (cid:80) i,j ∈ I p i ˆ P t,τ,i,j (cid:80) i ∈ I p i where ˆ P t,τ,i,j = p j P t,τ,j,i /p i then in [5] it is shown that ρ ( A ) = ˆ ρ ( A ). A major consequenceof this relation is that the total cost function ρ ( A , . . . , A q ) remains unchanged under time MEAGAN CARNEY*, HOLGER KANTZ † reversal and without loss of generality we may replace P t,τ,i,j with the time reversible matrixwith entries, R t,τ,i,j = 12 (cid:18) P t,τ,i,j + p j P t,τ,j,i p i (cid:19) (5)The authors in [6] have shown that solving a relaxation of the min-cut problem producesnear optimal values of the total cost function. Here, the greedy algorithm looks for apredefined number of q partitions on the graph defined by the Laplacian corresponding tothe symmetric matrix ( P t,τ,i,j + P t,τ,j,i ) /
2. On the other hand, it is shown in [5] that nearoptimal values can be obtained by performing a weighted fuzzy clustering algorithm onthe space generated by the eigenvectors of the time reversible matrix R t,τ directly.We take a different approach by performing spectral clustering equipped with K -meansto the reversible matrix R t,τ . This approach combines the benefits of both methods in-troduced in [5] and [6]. The Laplacian of the matrix R t,τ ensures the eigenvalues arepositive and real. The eigenvalues near zero and corresponding eigenvectors form thelower-dimensional subspace upon which the vectors of R t,τ are projected. This projectioncan be understood as a lower-dimensional representation of R t,τ without much loss of in-formation. On the other hand, choosing the matrix R t,τ over the symmetrized P t,τ matrixprovides spectral information on the number of expected almost-invariant sets and hencethe appropriate K value prior to clustering.We use the unnormalized graph Laplacian of the time reversible matrix R t,τ defined by, L t,τ = D t,τ − R t,τ where D t,τ is the diagonal matrix with entries D t,τ,i,i = (cid:80) i R t,τ,i,j . The matrix L t,τ issymmetric and positive semi-definite with 0 = λ ≤ λ ≤ · · · ≤ λ m . Moreover, solving forthe optimal q partition of L t,τ is equivalent to solving a relaxation of the min-cut problem[17].The spectral approach to partitioning L t,τ involves projecting the eigenvectors of L t,τ onto a smaller (cid:96) -dimensional subspace made of the ”most important” eigenvectors of L t,τ .Provided the perturbation ε is sufficiently small to produce a stochastic, irreducible matrix R t,τ , L t,τ has one eigenvalue equal to 0 and q − q (almost) invariant sets under the flow f ε [6, Thrm. 7.1]. Hence, the subspace sp { v , . . . , v q } can be used to identify the q almost invariant sets of P t,τ where (cid:96) = q − q − L with nearly optimalresults from clustering algorithms [9]. In our application, we choose to project onto sp { v } ,that is an (cid:96) = 1-dimensional subspace made of v corresponding to the second smallesteigenvalue of L t,τ .Given m blocks B i , each B i is represented by a point in R m +1 where the coordinate valueis given by the corresponding probability transition values to each block B j provided in R t,τ .The goal of K -means is then to find hyperplane separations between the clusters in R m +1 by minimizing the objective function defined by the sum of the inner cluster distances. Werefer the reader to Appendix A for a more rigorous explanation. Introducing L t,τ allowsus to run K -means using the projected, lower-dimensional representation of the vectors in R t,τ which often results in a set of more easily separable clusters. The result from K -means OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 7 clustering with K = q initial centroids will provide our optimal q spatial partition resultfor R t,τ . For a flow f ε representing movement in the atmosphere the partitioning result for R t,τ provides the geographic location of almost invariant sets (or eddies) where movementis trapped over the time interval [ t, t + τ ] for τ on the order of a day.4. Importance Sampling : Finding Trajectories Likely to End in anAlmost Invariant Set
We will refer to the subset A = ∪ i ∈ I B i made up of B i ⊂ X = R as the almost invariantset over the time interval [ t, t + τ ] for τ on the order of a day. We will require that A beconnected and define the center of A as γ = (( x (min) − x (max)) / , ( y (min) − y (max)) / φ ( x ) = − log( d ( x, γ )) x ∈ X where d is the Euclidean metric so that φ ( x ) is maximized as it approaches the centerof the almost invariant set. Let X N = φ ◦ f ε ( x N, , t ; τ ), where x N, is the N th startingposition, be a sequence of random variables representing the value of our observable φ asa function of the end position of the N trajectories on X run under the flow f ε from t upto time t + τ . Remark.
The importance sampling methods outlined in this paper are suggested as a start-ing point for almost invariant sets of any topology; however, they work best for symmetric,convex sets by the definition of the center γ . It is possible to obtain varying results fordifferent choices of γ even for nonsymmetric convex sets; for example, consider the set A as an ellipse and γ as either one of the focus points. Our justification for using this methodrelies on the fact that GPA results in a spatial distribution of points only with the highestdensity near γ . Remark.
Allowing φ ( x ) = − log d H ( x, S ) where S = { x : x ∈ δA } and d H = min d ( x, S ) isthe Hausdorff distance to S would be an interesting and more physically relevant choice formany climate questions. We have tested this definition and have found that it introducesseveral new complexities which future work would have to overcome including: (1) thedistribution of X ( x ) = φ ◦ f t ( x ) is numerically unclear with little theoretical support and(2) GPA would inevitably force surviving trajectories to group near the boundary.Genealogical particle analysis (GPA) is an importance sampling method that uses weightsto perform a change of measure in a reversible way so that rare events are sampled moreoften. These weights can be thought of as measuring the performance of the trajectory atspecified sampling times. Large values of the weight function imply that the trajectory isbehaving as if it comes from the target distribution. These trajectories will be cloned whilelow weight values indicate a trajectory that will be killed. Importance sampling algorithmsare often used to lower relative error of tail probability estimation because the change ofmeasure provides a set of trajectories that are more likely to end in a rare event. In ourcontext, running GPA will provide a pool of trajectories that are most likely to end in ouralmost invariant set over the time interval [ t, t + τ ]. MEAGAN CARNEY*, HOLGER KANTZ † One difficulty with GPA is determining a weight function that will change the measurein an appropriate way so that rare events are sampled more often. This choice depends onthe distribution of X N . We describe the explicit method used in this paper, however, werefer to [3, 15, 18] for more details on the general case. We assume that X N is distributedaccording to a unimodal distribution with tails decaying fast enough ( o (1 /N δ )) to zero.GPA has been successfully implemented numerically for random variables with symmetric,heavy-tailed distributions [3, 15, 18]. This assumption is supported by the numericalbehavior of the distribution of X N (see Figure 1). Remark.
Preliminary work using GPA on the Lorenz ’63 model equipped with the observ-able φ ( x ) = − log( d ( x, γ )) where γ ∈ X is some fixed point in the space also suggests X ( x ) = φ ◦ f t ( x ) follows a unimodal distribution with fast-decaying tails.1. Initiate N = 1 , . . . , M starting trajectories uniformly distributed over the space X .2. For n = 1 . . . (cid:98) τ /T (cid:99) where τ is the final integration time. T is referred to as the resampling time . Remark.
It is important to balance T between the correlation time of X N andthe Lyapunov time. Values of T taken too small can result in highly correlatedtrajectories (many of clones of a single trajectory) while too large can result in arelaxation back to the original distribution.2a. Iterate each trajectory from time t n − = T ( n −
1) to t n = T n .2b. At time t n , stop the simulation and assign a weight to each trajectory x n,N given by, W N,n = exp( C ( φ ◦ f ε ( x n − ,N , t n − , t n ) − φ ◦ f ε ( x n − ,N , t n − , t n − ))) Z n where Z n = 1 N M (cid:88) N =1 W N,n c N,n = (cid:98) W N,n + u N (cid:99) where (cid:98)·(cid:99) is the integer portion and u N are random variables generated froma uniform distribution on [0 , M n = M (cid:88) N =1 c N,n
Clones are used as inputs into the next iteration of the algorithm. For large N,the normalizing factor ensures the number of particles N n remains constant;however, in practice the number of particles fluctuates slightly on each iteration n . To ensure N n remains constant it is common to compute the difference∆ N n = N n − M . If ∆ N i >
0, then ∆ N i trajectories are randomly selected OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 9 (without replacement) and killed. If ∆ N i <
0, then ∆ N i trajectories arerandomly selected (with replacement) and cloned.3. The final set of trajectories ˜ X N = φ ◦ f (cid:15) ( x τ/T − ,N , t τ/T − , τ ) tends to a new distri-bution as N → ∞ exponentially tilted by the constant C . (x) p r obab ili t y C = 0C = 1C = 2C = 2.5
Figure 1.
Exponentially tilted distributions for different C values underGPA. Densities are estimated with a normal kernel. C = 0 corresponds tothe original distribution.The probability of sampling events for increasing values of φ ◦ f ε ( x, t ; τ ) under the ex-ponentially tilted distribution increases for increasing values of C . Since φ is maximizedat γ , the set of end trajectories is the set that has the highest probability of entering orremaining in the almost invariant set over the time interval [ t, t + τ ]. For a flow f ε rep-resenting movement in the atmosphere this set of trajectories form the paths of particlemovements which are most likely to enter or remain near the center of a target eddy.5. An Application to the Portable University Model of the Atmospherefor Pollution Movement
Portable University Model of the Atmosphere (PUMA) and Reconstructingthe Atmospheric Flow.
The PUMA module is a subset of PlaSim, a planet simulationmodel of intermediate complexity provided by the Universitat Hamburg MeteorologicalInstitute. Like most atmospheric models, PUMA is a simplified model derived from theNavier Stokes equation in a rotating frame of reference. In contrast to a full atmosphericgeneral circulation model (GCM), moist processes such as humidity, cloud water, evap-oration and the like are omitted. For more information on this model we refer to [13].Among the many available output options in PUMA, one is the time-dependent northern † and eastern wind velocities for a user-defined set of latitude and longitude pairs taken ata number of different atmospheric levels.Our goal for this section is to reconstruct a 2-dimensional time dependent atmosphericflow model. We begin with a set of time series of northern and eastern wind velocities eachtaken at a single (lat, lon) pair (64 latitude and 128 longitude in total) making up our wholespace Y and 1-day intervals with geopotential height equal to 1,000 hPA (approximately100 meters above sea level). We will occasionally refer to a set of (lat, lon) pairs as the grid and a single (lat, lon) pair as a grid point . Values of the wind velocity are given initiallyin m/s. We convert this value to km/day. We create a new set of time series with valuestaken at 0.1-day intervals by linearly interpolating the km/day time series at each gridpoint over time. Let V ( (cid:126)x ( t ) , t ) = (cid:20) V x ( (cid:126)x ( t ) , t ) V y ( (cid:126)x ( t ) , t ) (cid:21) where V x ( (cid:126)x ( t ) , t ) is the eastern velocity magnitude and V y ( (cid:126)x ( t ) , t ) is the northern velocitymagnitude for (cid:126)x ( t ) = ( x ( t ) , y ( t )) points on the grid at fixed time t . Then a trajectory (cid:126)x ( t ) := F t ( (cid:126)x , t ; t ) is a solution to the non-autonomous ODE d(cid:126)xdt = V ( (cid:126)x ( t ) , t ) and F t ( (cid:126)x, t, τ )is the end position of a trajectory in the 2-dimensional space of latitude longitude pairsbeginning at (cid:126)x , time t and flowing for τ time. F t can be understood as some ε perturbation f ε of the deterministic flow f in the previous discussion.We numerically approximate F t by the following Runge-Kutta method. We first notethat, by design, there is a fixed vector field at every 1 / th of a day. Hence, we mayintegrate over the fixed vector field at each discrete time step t n where t n = 0 . n day(s)from t (starting time of the flow), T = 0 .
1, and n = 1 , . . . , (cid:98) τ /T (cid:99) . For this application, wechoose t = 360 and τ = 4 days.PUMA and many standard GCMs require time for the system to reach an state wherethe initial conditions of the model (temperature gradients, CO concentrations, radiation,etc.) provide meaningful and accurate time-series outputs (velocity fields, temperature,etc.). In a dynamical sense, this is time it takes for a trajectory to move from its initialposition to the (chaotic) attractor. We choose a starting time t ≥
360 days to ensurethe model has reached this state. Furthermore, we choose our specific t based on a pre-determined guarantee that (at least) one eddy will pass through our region of interest byrunning a simulation of the velocity field movement over several different time-intervals.Let F t ( (cid:126)x n − , t n − , t n ) = (cid:126)x n = (cid:20) x n y n (cid:21) represent the discretized flow and calculate, (cid:96) = V ( (cid:126)x n , t n ) = (cid:20) V x ( (cid:126)x n , t n ) V y ( (cid:126)x n , t n ) (cid:21) = (cid:20) V x ( x n , y n , t n ) V y ( x n , y n , t n ) (cid:21) (cid:96) = V ( x n + 0 . (cid:96) (1) , y n + 0 . (cid:96) (2) , t n ) (cid:96) = V ( x n + 0 . (cid:96) (1) , y n + 0 . (cid:96) (2) , t n ) (cid:96) = V ( x n + 0 . (cid:96) (1) , y n + 0 . (cid:96) (2) , t n ) OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 11
And, (cid:126)x n +1 = (cid:126)x n + 0 . (cid:96) + 2 (cid:96) + 2 (cid:96) + (cid:96) ) / (cid:126)x n +1 is the end position of the trajectory starting at (cid:126)x n and flowing from time t n to t n +1 . For those end trajectories (cid:126)x n that lie between the original grid points ( x, y ), welinearly interpolate the vector field V over the ( x, y ) grid at the fixed time t n to obtain V ( (cid:126)x n , t n ). Computing the Transition Probability Matrix over a Fixed Time Interval.
Webegin by converting the latitude longitude coordinates to kilometers from the origin (equa-tor, center) so that integration over the km/day velocity is appropriate. We create a finergrid over the box 36 ◦ N − ◦ N and 11 ◦ W − ◦ E corresponding to Europe at ◦ intervals.We will refer to this box as X ⊂ Y and each box made of four corners of (lat, lon) pairsinside X as B i i = 1 . . . m . At this resolution, m = 4 ,
896 making up a 68 ×
72 size grid ofboxes. In each B i , we take k = 25 uniformly distributed points (cid:126)b ( k, i ) with coordinates (inkm) at (cid:126)b ( k, i ) = ( x ( k, i ) , y ( k, i )) and integrate to find the path over a fixed time interval[ t, t + τ ] using the Runge-Kutta method described above. Remark.
We are interested in restricting our search for almost-invariant sets to the smallregion of the atmosphere over Europe and investigating those sets that may be the result ofmultiple eddy interactions. To locate individual eddies in an expanded region, it is often ofinterest to first use the Okubo-Weiss criterion as a Eulerian way of quickly approximatingthe location of an eddy based on a combination of the vorticity, stretching, and shearingrate of the velocity field. After the location of an eddy is estimated, one restricts numericsto the identified region and uses the transfer operator approach to find the almost invariantset based on the Lagrangian dynamics inside the restricted region [14].We compute the transition probability P t,τ,i,j as the number of points (cid:126)b t ( k, i ) that start(at time t ) in B i and end (at time t + τ ) in B j [6]. We refer to Figure 2 for an illustrationof this flow. Note that the (cid:80) mj =1 P t,τ,i,j is not necessarily equal to 1 because it is possiblefor (cid:126)b τ ( k, i ) to land outside of the large box X under the flow. Following the work of [14] weintroduce a mixing ”tile” into this system defined as a box B m +1 in which all b τ ( i, k ) landingoutside of X are injected with equal probability into any B i in the next step. Addition ofthis tile is done at the end of computation of the transition probability matrix P t,τ = P t,τ,i,j by adding a column at the i = m + 1 position with value(s) equal to 1 − (cid:80) mj =1 P t,τ,i,j foreach j and adding a row at the j = m + 1 position with value(s) equal to 1 / ( m + 1) foreach i .Next, the probability measure on X = ∪ mi =1 B i is approximated by (3) where, p i = m ( B i ) m ( X ) = Area of B i Area of X so that the cost function is exactly as stated in (4), that is for any A = ∪ i ∈ I B i , I = { , . . . , m } ρ ( A ) ≈ (cid:80) i,j ∈ I p i P t,τ,i,j (cid:80) i ∈ I p i † Hence, we may solve for the almost invariant sets of the system by looking at the corre-sponding time-reversible matrix R t,τ given by (5) with entries, R t,τ,i,j = 12 (cid:18) P t,τ,i,j + p j P t,τ,j,i p i (cid:19) Figure 2.
Starting boxes B i (black) and particles under the flow F t at time t + τ . Velocity field at time t + τ is indicated by arrows. Spectral Clustering of the Time-Reversible Matrix.
We perform spectral clusteringof R t,τ to find the almost invariant sets in the box X defined as the European spatial regionin the following steps.1. Form the unnormalized Laplacian of the matrix R t,τ , L t,τ = D t,τ − R t,τ where D t,τ is the diagonal matrix with entries D t,τ,i,i = (cid:80) mi =1 R t,τ,i,j .2. Choose the first (cid:96) eigenvalues λ , . . . , λ (cid:96) and corresponding eigenvectors v , . . . , v (cid:96) of L t,τ .3. Form a subspace made of S (cid:96) := sp { v , . . . , v (cid:96) } and project the m -dimensional rowvectors of L t,τ onto S (cid:96) .4. Run K -means on the projected (cid:96) -dimensional row vectors with a predetermined K value. K essentially tells the algorithm how many almost invariant sets areexpected. OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 13
Remark.
Although the first q eigenvalues λ , . . . λ q ≈ L t,τ can givea rough estimate on the number of almost invariant sets in the space (and hence K ), this may result in finer level sets of the same ”parent” almost invariant set.The boxes B i which form the partition of the space X are now represented by m -dimensionalpoints in R m . The i th row is an m -dimensional vector with entries corresponding to thetransition probabilities between B i and every other box in the space. Ideally, the spectralmethod finds groups of boxes B i which are most likely to have transitions within groupsand least likely to have transitions between groups. Physically, this is equivalent to findingregions in the atmosphere where trajectories move within the region but are unlikely tomove outside of the region. This result may fail if K > q the number of required regionsis forced to be larger than the number of true almost invariant sets which is why it issuggested that the eigenvalues of P t,τ , R t,τ , and L t,τ are considered before choosing thisvalue.Over the subspace X = 36 ◦ N − ◦ N and 11 ◦ W − ◦ E representing a portion of theatmosphere over Europe (1,000 meters above sea level), we find two connected almostinvariant sets indicated by a single cluster in the spectral clustering outcome. We choosethe set corresponding to Europe for importance sampling. See Figure 3 for an illustrationof these sets. Figure 3.
Two atmospheric eddies as almost invariant sets under the flow F over [ t, t + τ ]. Regions are indicated by different colors. Center of thechosen target eddy is marked with a red X. † Importance Sampling and Finding Trajectories Likely to End in an Atmo-spheric Eddy.
We initialize a set (cid:126)x ,N of M = 12 ,
000 uniformly distributed points overthe larger domain 30 ◦ N − ◦ N and 20 ◦ W − ◦ E . We will refer to this domain as Z andemphasize that X ⊂ Z ⊂ Y . We refer to Figure 4 for an illustration of these regions. Theobjective of importance sampling is to enrich the data set of trajectories which are likelyto end near the center γ ( A ) of the almost invariant set A . By definition, this set has somemixing with the external space. This is illustrated in Figure 5. Figure 4.
Grid of B i boxes where X = ∪ mi =1 B i shown in blue. Initialvalues for importance sampling are taken uniformly over the region Z shownin black.For ease of computation, we choose the resampling time T = 0 . t n = T n for n = 1 , . . . , τ /T is equivalent to the step-size of the Runge-Kutta approximation of F t ;however we may theoretically choose any value T ≥ . φ ◦ F t ( (cid:126)x ( t ) , t ) is near zero at t = 0 . T . For pragmatic reasons, weadd a constant value D to the observable so that φ ( x ) = − log( d ( x, γ )) + D ≥
0. This shiftby D ensures that negative values in the exponent of the weight function are the result ofa true decrease of the observable value from the previous step.At each step n of the importance sampling algorithm outlined in Section 4, the weightfunction W N,n = exp( C ( φ ◦ (cid:126)x N,n − φ ◦ (cid:126)x N,n − )) Z n OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 15
A Z/AA AZ/A AZ/A Z/A
Figure 5.
End position of trajectories under the flow different colors illus-trate regions of mixing: (green) starting from the almost invariant set A andending in the external space Z/A , (purple) starting in A and ending in A ,(red) starting in Z/A and ending in A and (blue) starting in Z/A and endingin
Z/A . The center of the almost invariant set is marked with an X.is applied to exponentially tilt the original distribution by C where F t ( (cid:126)x N,n − , t n − , t n ) = (cid:126)x n is the end position of the N th trajectory under the (numerically approximated) flow F t beginning at (cid:126)x N,n − , time t n − and running until time t n . We choose the tilting value C = 2 . γ (see Figure1 for an illustration) while higher values of C force almost all trajectories toward the tailand result in a break-down of the distribution. This phenomenon in importance samplingmethods is studied in numerical detail in [3]. The algorithm is run over the interval[ t, t + τ ] = [360 , location distribution of surviving trajectories after genealogicalparticle analysis (GPA). Of course, most surviving trajectories are located near the center ofthe almost invariant set. We note that surviving trajectories with clearly low probabilitiesof entering the almost invariant set are due to the distributional aspect of GPA. That is,there is still a probability, though much lower, of obtaining trajectories far from the center. † Remark.
Following previous remarks from Section 4, we observe that trajectories thatbegin near the center γ ( A ) of the almost invariant set can be sent further away from γ ( A )under the flow but still remain inside the almost invariant set. These trajectories are oftenkilled in importance sampling because the value of the observable φ = − log( d ( x, γ ( A ))) istoo low toward the edges of the almost-invariant set. This is illustrated in Figure 7 (a).Decreasing the tilting value C decreases the number of killed trajectories that start in thealmost invariant set; however, this adjustment also decreases (at possibly a different rate)the number of killed trajectories in the general space. In this investigation, we focus ontrajectories that end in a reasonably large ball about γ ( A ); however this phenomenon isour motivation for considering some modification of the Hausdorff distance in future work.Since ˜ X N = φ ◦ (cid:126)x N,n is maximized at the center of our invariant set, the set of backwardreconstructed trajectories provides a collection of initial points that are most likely to enterand/or remain in the almost invariant set A over the time interval [ t, t + τ ] = [360 , Figure 6.
End position of trajectories (blue) under the flow and (red) aftergenealogical particle analysis. Note that most surviving trajectories end nearthe center of the almost invariant set (marked with an X).
OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 17
We may use this set of trajectories to measure the probability of the region E sending aparticle near the center γ of A under the flow by defining, p E ( E → B ( γ ( A ) , r )) = { x ,n ∈ E : x N,n ∈ B ( γ ( A ) , r ) } { x N,n ∈ B ( γ ( A ) , r ) } where we take E ⊆ Z and B ( γ ( A ) , r ) ⊆ Z is the ball centered at γ ( A ) with radius r . In thisway, we can determine the region that has the highest probability of trajectories ending upin A . We can also determine the most probable paths and set of initial conditions a regionwill take to get to A . Figure 8 shows the set of 10 ◦ × ◦ boxed regions E color coded interms of this probability. Application of this method can be used to determine the regionswith the highest probability of pollution exchange with a target eddy.6. An Application to the Portable University Model of the Atmospherefor Storm Tracking
We have shown in the previous example that the methods outlined in this paper can beused to track the collection of particles (or pollution) in the atmosphere to a fixed almostinvariant set that is formed by the background movement of multiple eddies through thespace. Now we consider a time-dependent almost-invariant set formed from a single eddy.Storm systems, such as hurricanes, have Lagrangian dynamics similar to that of an almost-invariant set so we can use these naturally occurring atmospheric eddies as a foundationalmodel for storm movement.We investigate a time-dependent almost invariant set as a function of a shorter timestep and use importance sampling to find its most likely path. Each step in the path isdetermined by a transition probability matrix, built over a small time interval on whichthe almost invariant set is defined, with states given by a spacial grid. Transitions aretaken as the probability of trajectories starting in a region and ending near the center ofthe almost-invariant set. Estimates of these probabilities are found by using genealogicalparticle analysis to enrich the set of trajectories likely to end near the almost-invariant set.As a rule, all notation in this section is carried over from the previous example.We numerically approximate the flow F t , built from the same northern and easternvelocity field outputs of PUMA, by the Runge-Kutta method described in Section 5. Next,we use the method described previously to find the almost invariant set in the region X ⊂ Y over J = 3 non-overlapping, consecutive time windows of length equal to 1 day,[ t + J, t + J + τ ] = [360 + J,
360 + J + 1]. The result is a time-dependent almost invariantset found over three discrete time intervals; one set is found for each time interval. Thelength of the chosen time intervals is relative to the movement speed of the almost invariantset. Time windows of a shorter length do not show a significant amount of movement ofthe almost invariant set while time windows of a longer length produce overlapping eddiesresulting in almost invariant sets of a different form. Our result is an almost invariant setthat moves spatially as a function of discrete time without much change in geometry.For each J , we run genealogical particle analysis on the set of uniformly distributedparticles over Z ⊂ Y with starting time t + J and termination time t + J +1. Recall that GPAreturns a set of trajectories that behave as though they come from the exponentially tilted † distribution where there is a higher likelihood of obtaining larger values of the observable φ ( x ) = − log d ( x, γ ( A ( J ))) where γ ( A ( J )) is the midpoint (center) of the J th correspondingalmost invariant set. Hence, the outcome is the set of trajectories most likely to end near γ ( A ( J )). The resampling time is again taken at T = 0 . t n = T n , n = 1 , . . . , (cid:98) τ /T (cid:99) = 10. Backwards reconstruction of surviving trajectories is then used todetermine the set of initial points which are most likely to end near γ ( A ( J )).In this example, the set of possible starting regions is defined after GPA as the setof 5 ◦ × ◦ boxes covering all of Z . For each J , we have an associated region E J ⊂ Z corresponding to the starting region that has the highest proportion of initial points fromsurviving trajectories (over all starting regions). Since E J has the highest probability ofsending trajectories near γ ( A ( J )) at time t + J + 1, this region should provide us with themovement direction of the almost invariant set A ( J ) → A ( J + 1) (and its correspondingcenter γ ( A ( J )) → γ ( A ( J + 1))) defined over [ t + J + 1 , t + J + 2]. Using each of the J = 3invariant sets found previously from the PUMA flow approximation, we illustrate in Figure10 that E J can provide some reasonable indication of movement direction for the almostinvariant set in the next time step. Remark.
The user-chosen resolution of the grid of regions E depends on the sparsity of theGPA outcome and the original resolution of the velocity field. This choice of resolutioncan produce slightly varying results for the predicted direction of movement.7. Conclusion
Almost invariant sets in the atmosphere are physically interesting because they relateto single eddies and eddy interactions. In this investigation, we look at almost invariantregions of atmospheric flow represented by the numerically approximated flow taken fromvelocity fields of the Portable University Model of the Atmosphere. We use a modifiedset of tools revolving around a well-studied transfer operator approach to estimate theseregions. In particular, we approximate the Perron-Frobenius operator by the transitionprobability matrix for the flow over the European subregion and use spectral K -meansclustering to find almost invariant sets only located over Europe.It can be seen for longer time intervals that the almost invariant sets are formed bymultiple eddy interactions such as particles trapped between the spin of two eddies; whereasshorter time intervals have almost invariant sets corresponding to a single eddy. For theformer, one can ask questions about the regional origin of trajectories ending inside thealmost invariant set and their most probable paths. For the latter, one can ask questionsabout the path of such an eddy in the space. To study these trajectories, we employ awell-known importance sampling algorithm, called genealogical particle analysis, not usedin this context to-date.Current literature has focused on using importance sampling methods to decrease therelative error of an estimated rare event probability by forcing rare events to occur morefrequently. We show that these methods can also provide useful information on the setof trajectories likely to end in an extreme event. For the interest of this study, we haveintroduced an observable that defines the extreme event as being near the center of an OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 19 almost invariant set. In this setting, we show that the surviving trajectories obtainedfrom importance sampling can provide information on probable paths and initial regions oftrajectories that end in an almost invariant set under an atmospheric flow. We complete ourinvestigation by motivating and illustrating some important examples where informationabout trajectory movement toward the center of an almost invariant set in the atmosphereis useful and physically relevant: storm movement and origin of pollution.In future work we plan to apply these techniques to real hurricane data where the fixedpoint γ may be taken as some point outside of the time-dependent almost-invariant set.The outcome of importance sampling would then give us the probability of a hurricanemoving over a given region. It would also be interesting to consider importance samplingmethods for the sequence of maxima M n = max { X , . . . , X N } where X N = φ ◦ f ( x N ).This would limit the set of original distributions to the family of generalized extreme valuefunctions and possibly provide a new way of using the Hausdorff distance in the definitionof φ ( x ) = − log( d ( x, γ )). Furthermore, a complete shift of the generalized extreme valuedistribution under exponential tilting would result in a higher density around γ and lessuniformly distributed points about the whole space. † Figure 7. (a) Parent initial positions of surviving trajectories after ge-nealogical particle analysis. Note the cluster in the southwest region of thespace. This distribution is used to determine the regions with the high-est likelihood of ending near the center γ ( A ) of the almost invariant set A (marked with an X). (b) Initial positions of surviving trajectories ending insome ball B ( γ ( A ) , r ) of radius 10 ◦ and not beginning in the almost invariantset. (c) Uniformly sampled initial positions (red) of (b) with full trajectories(blue). OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 21
Figure 8.
Surviving trajectories after genealogical particle analysis (blue).The set of all surviving trajectories x N,n ∈ B ( γ ( A ) , r ) and x N, / ∈ A with r = 10 ◦ are marked in red and their parent initial trajectories x ,n (black).Grid boxes are 10 ◦ × ◦ with darker orange indicating a higher p ( E → B ( γ ( A ) , r ) value. The regions with the highest probability of ending near γ ( A ) are 10 ◦ E − ◦ W and 30 ◦ N − ◦ N . † Figure 9.
Almost-invariant sets found from the transfer operator methodof J = 3 non-overlapping, consecutive time windows of length 1 day. Timeintervals [ t + J, t + J + τ ] = (a) [360 , , , OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 23
Figure 10.
Movement prediction of the almost invariant set. The centerfor step J is marked with a red X and the corresponding initial positionsof the end surviving trajectories within B ( γ ( A ( J )) , r ) are indicated by redpoints. The J +1 almost invariant set is highlighted in gray, its correspondingvelocity field is represented by black quivers and center is marked with a blackX. (a) J = 1 and (b) J = 2. The movement direction probability is takenas the proportion of (surviving) initial positions in the regions marked bythe grid. These probabilities are highlighted in orange with darker valuesindicating a higher probability. † Appendix A. K -Means The standard K -means algorithm for a set of m nodes n represented by m vectors in R H is given by,Step 1 Start with K random partitions P j of the space R H .Step 2 Compute the centroids (means) of these partitions as C j = (cid:80) n ( (cid:96) ) ∈ P j n ( (cid:96) ) / card( P j )where C j ∈ R H .Step 3 Assign n ( (cid:96) ) to the partition P j with the minimum (squared) euclidean distancebetween n ( (cid:96) ) and C j .Step 4 Update the algorithm by recalculating the centroids (means) of P j .The algorithm continues by repeating steps 3 and 4 until the assignments no longer change.This is equivalent to finding the steady state of the objective function given by,min K (cid:88) j =1 m (cid:88) (cid:96) || n ( (cid:96) ) − C ( j ) || R H , the minimum sum of the (squared) euclidean distances between each node and its assignedcentroid. OURCES AND SINKS OF RARE TRAJECTORIES IDENTIFIED BY IMPORTANCE SAMPLING 25
References [1] Carney, M., Azencott, R., Nicol, M.:
Nonstationarity of summer temperature extremes in Texas , IntJ Climatol. Robust regional clustering and modeling of nonstationary summer temperatureextremes across Germany , Advances in Statistical Climatology Meteorology and Oceanography. Analysis and Simulation of Extremes and Rare Events in ComplexSystems arXiv: 2005.07573. preprint.[4] Dellnitz, M., Junge, O.:
On the Approximation of Complicated Dynamical Behavior , JSTOR. (2)491-515, 1999.[5] Froyland, G.: Statistically Optimal Almost-Invariant Sets , Physica D. (3-4) 205-219, 2005.[6] Froyland, G., Dellnitz, M.:
Detecting and Locating Near-Optimal Almost-Invariant Sets and Cycles ,SIAM J. Sci. Comp. (6) 1839-1863, 2003.[7] Froyland, G., Padberg-Gehle, K.: Almost-Invariant and Finite-Time Coherent Sets: Directionality,Duration, and Diffusion , Bahsoun W., Bose C., Froyland G. (eds) Ergodic Theory, Open Dynamics,and Coherent Structures. Springer Proceedings in Mathematics and Statistics, Springer. , 2014.[8] Froyland, G., Schwalb, M., Padberg, K., Dellnitz, M.: A transfer operator based numerical inves-tigation of coherent structures in three-dimensional Southern ocean circulation , Proceedings of theInternational Symposium on Nonlinear Theory and its Applications, 2008.[9] Hastie, T., Tibshirani, R., Friedman, J.:
The Elements of Statistical Learning: Data Mining, Inferenceand Prediction
Almost sure convergence of maxima for chaotic dynamical systems ,Stochastic Processes and their Applications , 126, 3145-3170, 2016.[11] Junge, O., Marsden, J., Mezic, I.: Uncertainty in the Dynamics of Conservative Maps , 2004 43rdIEEE Conference on Decision and Control (CDC) (IEEE Cat. No.04CH37601). Extremes and Recurrence in Dynamical Systems , Wiley 2016.[13] Lunkeit, F. Blessing, S., Fraedrich, K., Jansen, H., Kirk, E., Luksch, U., Sielmann, F.:
Planetsimulatorusers guide version 15.0 , Meteorological Institute of the University of Hamburg, Hamburg, 2007.[14] L¨unsmann, B., Kantz, H.:
An extended transfer operator approach to identify separatrices in openflows , Chaos. Computation of extreme heat waves in climate models using alarge deviation algorithm , Proc. Natl. Acad. Sci. USA, (115) 24-29, 2018.[16] Vallero, D.: Fundamentals of Air Pollution 5th ed. , Elsevier, 2014.[17] von Luxburg, U.:
A Tutorial on Spectral Clustering , Statistics and Computing. (4), 2007.[18] Wouters, J., Bouchet, F.: Rare event computation in deterministic chaotic systems using genealogicalparticle analysis , J. Phys. A.,49