A closing scheme for finding almost-invariant sets in open dynamical systems

We explore the concept of metastability or almost-invariance in open dynamical systems. In such systems, the loss of mass through a ``hole'' occurs in the presence of metastability. We extend existing techniques for finding almost-invariant sets in closed systems to open systems by introducing a closing operation that has a small impact on the system's metastability.


Introduction
The classical setting for dynamical systems is one in which a transformation T : X → X is iterated to generate forward trajectories of infinite length.Interesting phenomena occur if a "hole" H ⊂ X is introduced to X and any trajectory falling into H immediately terminates and is no longer considered.We will call a dynamical system open if its domain contains holes and closed if not.The notion of an open dynamical system is not new and we do not attempt to exhaustively cover all of the literature on open systems.Our focus is on systems with some (piecewise) smoothness and their ergodic properties.In recent years, there has been a burgeoning interest in the mathematical characterisation of these types of open systems [36,17,34,32,7]; see also the survey article [18] and references therein.Early work on the ergodic theory in one-dimensional dynamics was carried out by Pianigiani and Yorke [39], and further results have been established for a range of one-dimensional maps [9,10,11,6].Throughout we will assume that X is a smooth Riemannian manifold embedded in R d and T a piecewise smooth transformation on X.
We now begin by describing three simple mental models for open, metastable, and openmetastable systems.In these simple models, X is "coarse-grained" (partitioned) into two or three large subsets.
Example 1.1.Coarse-grain the state space X into two disjoint parts; the hole H and its complement, denoted by X r (r for "restricted").A schematic representation of the coarsegrained dynamics is shown in Figure 1.Given some natural reference probability measure, at this coarse-grained level, we can represent the probabilities of exchange between the two parts via a 2×2 stochastic transition matrix P open , where P open,ij is the conditional probability of moving from state i to state j.For open systems, the transition matrix would have the following form: where we think of ǫ X r →H > 0 as small.Applying the coarse-grained probabilistic dynamics given by P open , the amount of material remaining in X r after k steps is proportional to A second phenomenon of recent interest for closed systems is metastability and the characterisation of almost-invariant sets.Metastable systems have dynamical processes operating on two or more different timescales.Typically there are fast processes which mix the state space and produce fast decorrelation of observables, and there are slow processes that often control switching between different (larger) parts of the state-space.A simple schematic diagram at a coarse-grained level is shown in Figure 2; one imagines that within X 1 and X 2 there is fast mixing compared to the mixing between X 1 and X 2 .
Example 1.2.Coarse-grain the state space of a closed dynamical system into two disjoint parts, denoted X 1 and X 2 , where X 1 ∪ X 2 = X.The system is metastable if the conditional probabilities of exchange between the two parts are both small.The transition matrix representing the coarse-grained conditional transition probabilities has the following form: where ǫ X 1 →X 2 > 0 and ǫ X 2 →X 1 > 0 are both small.
The coarse-grained regions X 1 , X 2 are known as almost-invariant sets or metastable sets and are usually defined as regions of state space between which there is limited exchange or in which trajectories tend to stay for a long time.In smooth dynamics, the dynamical creation of almost-invariant sets may be due to co-dimension 1 stable and unstable manifolds, for example, as in lobe dynamics [41,40,47], and there is numerical evidence that there is a relationship between the geometry of invariant manifolds of low period orbits and dominant almost-invariant sets [29].However, in general, determining the location of highly almostinvariant regions is a difficult problem.Transfer operator approaches have arguably been the most successful in applications: identifying molecular conformations as almost-invariant sets [43], locating oceanic gyres as almost-invariant sets [30,15], studying noisy epidemic models [3], identifying "ghost-rods" in fluid mixing via almost-invariant sets [44], and in the time-dependent setting, determining the three-dimensional location of the Antarctic polar vortex [31] as a coherent set.
There is a clear similarity between open and metastable systems: both systems give rise to quasi-stable behaviour on a subset of the state space before trajectories eventually leave that subset.The similarity between the two characterisations has been exploited in [34,7], both of which treated metastable systems as perturbations of nonergodic systems, and in [32], which uses an approach, originally created to identify almost-invariant sets, to partition closed systems into two open systems with slow escape rates.The distinguishing feature of open systems is that the trajectories cannot re-enter the subset X r once they have fallen through the hole, whereas in metastable systems, trajectories can exit and reenter X 1 or X 2 with some small probability.For open systems one refers to the escape rate from the restricted domain, whereas in metastable systems one refers to the exchange rate.The purpose of this paper is to study the metastability properties of dynamical systems in the presence of holes and specifically to develop techniques to identify and characterise almost-invariant sets in open systems.
Example 1.3.Coarse-grain the state space X into X 1 , X 2 and H.This coarse-graining produces an open metastable system if the transition matrix of conditional probabilities restricted to X r = X 1 ∪ X 2 has the form where all ǫ values are small and positive.A schematic diagram of the coarse-grained dynamics is shown in Figure 3.The matrix P open−metastable is sub-stochastic.
Computational approaches to determine almost-invariant and metastable sets in closed dynamical systems are typically based around numerically accessible discretisations of the Perron-Frobenius operator, eg.[20,13,16,19,14].These discretisations produce a transition matrix (similar to the coarse-grained matrices in the above examples, but with a much finer graining).The spectrum of this transition matrix is used to quantify the number of almost-invariant sets [19,14].The eigenvectors corresponding to large eigenvalues are used to heuristically locate almost-invariant sets.An improved heuristic for non-time-symmetric dynamical systems was proposed in [27]; an induced time-reversible Markov chain was constructed from the transition matrix and the maximality properties of the spectrum and eigenvectors of time-reversible Markov chains were exploited to obtain rigorous bounds on maximum almost-invariance ratios.Using the construction of [27] recent numerical studies in [29] have demonstrated strong connections between co-dimension 1 invariant manifolds of low period points and boundaries of almost-invariant sets.In this paper our main contribution is to extend the methodology of [27] by applying a "closing operation" to the open system that has a small effect on the metastability properties.We state here our main result linking the metastability of a set A ⊂ X r under our new induced closed system, denoted by Wν (A), with the metastability of a set A ⊂ X r under the original open system, denoted by W r ν (A).
Theorem 1.4.Denote the restriction of T to X r by T r , and let ν be a conditionally invariant probability measure for T r (see definition in Section 2).Under mild conditions on ν, the following relationship holds: where H 1 = X r ∩ (T r ) −1 H is the set of points in X r that will fall into H in one step.
Thus, for holes with small ν-measure, the difference between the invariance (or metastability) ratios of the original open and induced closed systems is also small.An outline of the paper is as follows.Section 2 defines almost-invariance for closed and open systems and introduces the notions of invariant and conditionally invariant probability measures.Section 3 introduces the Perron-Frobenius operator in closed and open systems and describes the numerical discretisation scheme for approximating these operators.Section 4 contains our main results.It begins by summarising the content of [27], then introduces a closing operator to convert an open system into a closed system.Formulae for the discretised Perron-Frobenius operators describing this induced closed system are developed and we prove relevant quantitative relationships between the original open system and the induced closed system, such as Theorem 1.4.An illustrative example will be discussed throughout the paper.Section 5 describes numerical experiments for a two-dimensional system.Most proofs are deferred to the Appendix.

Invariant and conditionally invariant probability measures and metastability
In order to discuss the mass or probability transported by T we endow X with the Borel σ-algebra B and create the measure space (X, B, ℓ) where ℓ is Lebesgue measure normalised so that ℓ(X) = 1.For B-measurable T (for example, piecewise continuous T ), we call the dynamical system (X, B, ℓ, T ) a closed dynamical system.
for every A ∈ B(X).
Example 2.2.Suppose that we define a closed, piecewise linear map T : [0, 1] by: The parameter α can be varied between 0 and 1 2 , and for every value of α, Lebesgue measure ℓ on the unit interval is a T -invariant probability measure.When α = 0, there are two invariant sets, A ′ = [0, 1/2] and A = [1/2, 1].When 0 < α ≤ 1/2, (T, ℓ) is exact and thus mixing, and there are no nontrivial invariant sets (see eg. [46], Remark 1, p40).To emphasise the relationship between this map and the usual doubling map, we call this the α-shifted doubling map.We will set α = 1/16.The graph of this map on the unit interval is shown in Figure 4.
Throughout we suppose that the underlying dynamical system (X, T ) has a natural invariant probability measure µ.Under coarse-graining with respect to µ, one has ǫ X 1 →X 2 = µ(X 1 ∩ T −1 X 2 )/µ(X 1 ) and ǫ X 2 →X 1 = µ(X 2 ∩ T −1 X 1 )/µ(X 2 ).This motivates the definition of the invariance ratio: Definition 2.3.For a closed dynamical system with a T -invariant measure µ, let the invariance ratio for a subset A ∈ X be given by For metastable systems, we expect to find a coarse-graining By determining the X k that have large values we learn about the finite-time transport properties of T .
Let (X r , T r ) define an open dynamical system where T r := T | X r .The natural analogue of an invariant probability measure for T r is a conditionally invariant probability measure.
, and there is an obvious analogy with (2.1).We will later introduce a hole to the system of Example 2.2; numerically there appears to be a natural choice of conditional invariant measure ν, however there is no explicit analytic description of ν.We may similarly consider invariance ratios for subsets of X r under the action of T r .Definition 2.5.An open dynamical system with a T r -conditionally invariant probability measure ν, has an invariance ratio given by We have 1 the only mass that leaves A in one iteration immediately falls in the hole; there is no transmission to X r \ A. We will call a set 3 The Perron-Frobenius operator and discretising the dynamics We will assume that the invariant and conditionally invariant probability measures µ and ν are absolutely continuous1 with respect to ℓ and ℓ r , respectively, where ℓ r is normalised Lebesgue measure on X r .Then the Radon-Nikodym derivatives f := dµ/dℓ and f r := dν/dℓ r exist and are densities with respect to Lebesgue; ie.f, f r ≥ 0 and X f dℓ = 1, X r f r dℓ r = 1.

Perron-Frobenius operators
We will further assume throughout that T is nonsingular with respect to ℓ; that is, ℓ • T −1 is absolutely continuous with respect to ℓ.This allows us to define a linear operator describing the action of T on densities.
For the map T r , we assume One may also write P r (g) = P(½ X r • g).One has P r (f r ) = λ r 1 f r with λ r 1 > 0 and we call f r a conditionally invariant density.

Discrete Perron-Frobenius operators
As P and P r capture a significant amount of dynamical information we employ numerical methods to extract this information.We use a Galerkin-based procedure known as Ulam's method [45]; see also [37,26,4,22].We partition X = n i=1 B i into n disjoint connected sets of small diameter.In practical computations if X ⊂ R d the sets B i , i = 1, . . ., n are often d-dimensional cubes and henceforth we refer to them as boxes.We assume that the hole H is comprised of a finite union of boxes from this partition.Define the collection of all sets that are unions of boxes in X by B n , and similarly, the collection of all sets that are unions of boxes in X r by B r n .We also let I = {1, 2, . . ., n} , I r = {i ∈ I : For the open dynamics we similarly define π r : L 1 (ℓ r ) → span(½ B i : i ∈ I r ): recalling that ℓ r is normalised Lebesgue measure on X r , so ℓ r (A) = ℓ(A)/ℓ(X r ) for all B r n −measurable sets A. The discrete versions of P and P r are defined as: The matrix representations of Q and Q r (under multiplication on the left) acting on span{½ B i : i ∈ I} and span{½ B i : i ∈ I r } respectively are: and Note that Q and Q r are identical over i, j ∈ I r .Via a similarity transformation we define matrices P = M −1 QM and P r = (M r ) −1 Q r M r where M and M r are diagonal matrices with diagonal elements (M ii ) = ℓ(B i ), i ∈ I and (M r ii ) = ℓ(B i ), i ∈ I r , respectively.The matrix P is row stochastic and may be interpreted as a transition matrix; the entry is a conditional probability of a randomly chosen point in box i moving to box j in one iteration of the map T .Similarly, the matrix P r = (P r ) ij is defined with entries This matrix is row sub-stochastic and may be interpreted as a transition matrix for a transient Markov chain approximation of the open dynamical system T r .Using the leading left eigenvector p ≥ 0 (normalised so that i∈I p i = 1) of P , we will define an approximate absolutely continuous invariant probability measure (acim) for T by Similarly, an approximation of an absolutely continuous conditionally invariant measure (accim) ν is defined as where p r ≥ 0 is the leading left eigenvector of P r , normalised to sum to 1. Since P r is sub-stochastic, p r satisfies i p r i P r ij = λ r n,1 p r j , where λ r n,1 < 1 is the leading eigenvalue of P r .A rigorous justification for the approximation µ n goes back to Li [35] who considered expanding, piecewise C2 maps T of the unit interval.Since then, similar results, extending the class of maps have been developed [24,21,25,38].In the open setting, the approximation ν n has also been rigorously justified for expanding maps of the interval [1].
Example 3.4.Consider the α-shift map with H = [0, 0.0468].We partition the interval [0.0468, 1] into n = 9532 boxes each of width 10 −4 , compute 2 the transition matrix P r , and plot the left eigenvector, corresponding to the approximation of the accim, in Figure 5.Note that the invariant measure of the open system is quite different from Lebesgue.The leading eigenvalue of P r is approximately 0.9639.

Discrete invariance ratios for sets
Using the discrete measures µ n , ν n and transition matrices P, P r we now define discrete versions of W µ (A) and W r ν (A).Denote the index set of boxes that make up a set A ∈ B n by For A ∈ B r n we define a discrete version of W r ν (A) as We refer the reader to [28] for rigorous results concerning convergence of W µn (A) to W µ (A) as the box diameters go to zero.
Example 3.5.For the α-shift map, we compute the invariance ratios for intervals of the form [x, x + 0.5], where the left end point x is between [0, 0.5] for the closed system.This will enable us to compare the invariance ratios for all possible intervals of µ-measure 1/2.In practice, we begin by partitioning the unit interval into 10000 boxes.We compute the invariance ratio for the closed system with x = 0, increment x by 0.0001, and repeat the calculations, stopping when x = 0.5.The resulting invariance ratios are shown in Figure 6a.We find that the maximum value of W µn (A) for A in the class of intervals described above is W µn ([0, 1/2]) = W µn ([1/2, 1]) = 0.9375.To repeat this exercise for the open system, we discard the boxes corresponding to the hole H = [0, 0.0468], which leaves n = 9532.For comparability with the closed system, we take intervals of the form l(x) = [x, y(x)], where ν n (l(x)) = 1/2.We set x = 0.0468 and calculate y(0.0468),W r νn (l(0.0468)).We then increment x by 0.0001 and repeat, stopping when x = 0.6840 because ν n ([0.6840, 1]) = 1/2.The results are shown in Figure 6b.The maximum value of W r νn (l(x)) is W r νn (l(0.5001))= 0.6821.
For the α-shift example we find that the invariance ratio for the closed system is always greater than that for the open system for this particular class of intervals.This is at least in part because the interval for the open system contains the preimage of the hole H 1 = [0.7188,0.7421] whenever x > 0.1562.However, it is also because of the dynamics of T r and the difference between the measures µ and ν.To see why, consider the graph of the map, overlaid in Figure 6b, and as an illustrative example, consider the set l(0.0468).The set of points which will leave the interval (the 'exchange set') under iteration of T r are those between [0.5608, 0.6840], i.e., those whose image is outside of the original interval.But this exchange set has higher ν-measure than the rest of the interval (see Figure 5), and consequently the invariance ratio is lower.By contrast, the closed system has Lebesgue as its invariant measure, so it will never give additional weight to the exchange set.Note that this argument only holds because the denominators of the invariance ratios W µn (A) and W r νn (A) are the same, and we need only compare the numerators.
In Example 3.5 we considered only intervals of invariant measure 1/2.The next section describes a heuristic to find subsets of the unit interval of measure no greater than 1/2 with high invariance ratios.

An induced closing operation and time-reversal symmetry
We begin by briefly recapping the relevant constructions from [27] and then describe the construction of our induced closing and its properties.The quantity W µn (A) is time-symmetric even if the original dynamics given by T , or the discrete dynamics given by P are not.Mathematically, W µn (A) = i,j∈I A p i Pij / i∈I A p i , where Pij := p j P ji /p i is the transition matrix governing the Markov chain P in backward time (see e.g.[5] for details on P ).Thus, one can define R = (P + P )/2, a reversible Markov chain with invariant measure p, and obtain Utilising bounds due to [33], this construction yields rigorous upper and lower bounds on the invariance ratio of sets in B n : where λ R 2 denotes the second largest eigenvalue of R. The corresponding right eigenvector v of R is used in a heuristic procedure to construct sets A with high values of W µn (A).One chooses a "threshold value" c and defines sets A 1 (c) = i:v i ≤c B i and A 2 (c) = i:v i >c B i .The threshold c may then be varied to maximise W µn (A 1 (c)) or W µn (A 2 (c)) while maintaining µ n (A 1 (c)), µ n (A 2 (c)) ≤ 1/2.For further details, we refer the reader to [27] or [29].We wish to apply this methodology based on the reversible transition matrix R, as this approach has found success in physical oceanography [30,15], epidemic models [3], and 3D fluid mixers [44].
A direct application of the methodology in [27] to the open setting is not possible as there is no obvious candidate for Pr .We therefore will construct a "closed version" of P r which approximately maintains the metastability of properties of P r .We then derive its time-reversed counterpart.

The 'Closure' of the Open System
Define H 1 = X r ∩ T −1 (H), the set of points in X r that are about to fall into the hole.Further, define a functional τ : L 1 (ℓ r ) → R to be the expectation of g : X r → R on H 1 with respect to ℓ r : We will now construct new Markov operators for the original and the discretised spaces.
Definition 4.2.Define the operator P : L 1 (ℓ r ) → L 1 (ℓ r ) by We call P the closure of P r .The action of P is the same as that of P r , except that the part of the density g that is about to fall into H is redistributed over all of X r according to the conditional invariant density f r , scaled by the amount of mass about to fall into H.Lemma 4.3.P is a Markov operator with fixed point f r .Furthermore, if P r g = λg for g ∈ L 1 (ℓ r ) and λ = 1 then Pg = λg.
Proof.See appendix.
Thus, the eigenspectra and eigenfunctions of P r and P are identical, except for the leading eigenvalues.
We wish to define an Ulam-type approximation Q : span{½ B i : i ∈ I r } by Q := π r P. Define H 1 as the index set of boxes that intersect with the preimage of the hole, i.e.
Lemma 4.4.The matrix representation of Q on span{½ Proof.See Appendix.
The term p r j /ℓ r (B j ) may be interpreted as the discrete version of f r , while ℓ(B i ∩ H 1 ) is simply the ℓ-proportion of each box that will escape through the hole on the next iteration.Define P = (M r ) −1 QM r , and the vector τ n by τ n,i = 1 − j∈I r P r ij .
Lemma 4.5.The matrix P is row stochastic, has p r as a fixed left eigenvector and satisfies where τ n is a column vector with elements and p r is a row vector.
Proof.See Appendix.
We call P the closure of P r .The idea of resurrecting a Markov chain at time it leaves a particular set goes back to [2], but only recently has the idea been exploited fully (see for example [23,8]).In this case, since p r is the leading fixed left eigenvector, it follows that ν n , the accim for the open system, is an invariant measure for the discrete closure system.We would now like to define invariance ratios in terms of the closure operators P and P .It is not possible to do this in the same way that we defined invariance ratios for the closed and open systems, because those were defined in terms of the transformations T and T r , and there is no deterministic transformation that corresponds to the closure operators.However, the following lemma shows that we can write the invariance ratio W r ν (A) for A ⊂ X r directly in terms of the operator P r .For the remainder of the paper we assume that the conditional invariant density f r is strictly positive.Lemma 4.6.Let g 1 , g 2 ν denote the usual scalar product Proof.See Appendix.
Using this Lemma we can write W r ν (A) = ½ A , P r ν ½ A ν /ν(A).This permits a natural definition for the invariance ratios under the closure operators.
Definition 4.7.For a measurable set A ⊂ X r and Pν (g With this definition, the following theorem shows that the metastability properties of the system induced by the closure operation are closely related to those of the open system.Theorem 4.8.Suppose f r > 0. For measurable A ⊂ X r , one has Proof.See Appendix.
The above theorem shows that the invariance ratio for the closure of the open system is equal to the invariance ratio of the restricted system, plus an adjustment term of at most In analogy to (4.6), we define a discretised version of Wν (A) for sets A ∈ B r n : Definition 4.9.For A ⊂ B r n , we set In analogy to Theorem 4.8, one has Lemma 4.10.
Proof.See Appendix.
Example 4.11.For the 1/16-shifted doubling map with H = [0, 0.0468], we again calculate the invariance ratios for all intervals of the form l(x) = [x, y(x)], where ν n (l(x)) = 1/2 and the left end point x is between [0.0468, 0.6839].The interval [0.0468, 1] is partitioned into 9532 boxes of width 10 −4 , we set x = 0.0468, compute y(0.0468),W r νn (l(0.0468)) and Wνn (l(0.0468)),advance x by 0.0001, and repeat.Figure 7 shows that W r νn (A) is equal to Wν (A) when the set A does not overlap the preimage of the hole H 1 = [0.7188,0.7421]], cf Theorem 4.8.This is the case when x < 0.1562.Theorem 4.8 also dictates that the discrepancy between Wν (A) and W r νn (A) is bounded above by 1 − λ r n,1 = 0.0361.νn (A) for intervals of ν-measure 1/2 for the 1/16-shifted doubling map.The difference is zero when the interval does not include the preimage of the hole (which is demarcated with dashed lines), and is bounded from above by 1−λ r n,1 = 0.0361.

Utilising eigenvector maximality properties of R
Now that we have a stochastic matrix P (Lemma 4.5) and a relationship between W r νn (A) and Wν (A) (Lemma 4.10), we may apply the methodology of [27] to find highly almost-invariant sets for the induced closed system.
Let Pij = p r j Pji /p r i , i, j ∈ I r denote the transition matrix for the Markov chain P in backward time.The invariance ratio Wν (A) is time-symmetric: Lemma 4.12.For A ∈ B n , Proof.Straightforward to verify.
We define a time-symmetrised version of P as The Markov chain with transition matrix R is reversible and preserves the measure p r .Using the matrix R and the relationship between Wν (A) and W r νn (A) given in Lemma 4.10 we can produce rigorous bounds for almost-invariant ratios for the original open system.Theorem 4.13.One has the following bounds for the maximum invariance ratio in an open system: Proof.Apply (4.1) to the matrix R and invariant measure ν n .Theorem 4.8 yields the result.
Furthermore, we may apply the heuristic from [27] to the induced closed system to find sets with high invariance ratios for the original open system.
Example 4.14.To the α-shifted map with α = 1/16 we have introduced a hole H = [0, 0.0468].We compute P r with 9532 boxes.The leading two eigenvalues of P r are {0.9639,0.8920}.We form P ; its leading eigenvalue is 1 and the remaining outer spectra of P is the same as that of P r (cf.Lemma 4.3).
We then form R and compute its second eigenvalue as 0.9736 and corresponding right eigenvector v 2 ( R).We are also now able to compute the bounds in Theorem 4.13, and state that the maximal invariance ratio W r νn (A) must be in the interval [0.7343, 0.9868].We now wish to apply the thresholding algorithm described in Section 4, or in greater detail in Steps 6-7 of Algorithm 1 in Section 5 (see also [27] or [29]).Briefly summarising, we select a value c to maximise B i for some candidate eigenvector v. Our suggested approach is to take v to be the second right eigenvector of R.However, we will also compare this result with sets obtained from three other candidate methods that one might use in the absence of our new approach.Our goal in doing so is to see whether the new matrix R is an improvement over existing alternatives.
Method 1 Threshold the second right eigenvector of R (Figures 8a, 8b).   4 For non-invertible maps the right eigenvectors of P and P r are typically highly irregular and one instead uses left eigenvectors, which are smoother.This approach has been taken for closed systems in e.g.[16].8h).The values of W r νn (A opt ) and ν n (A opt ) resulting from Methods 1-4 are summarised in Table 1.It is clear that Method 1 is the superior approach.For comparison we may relate these results to those from Example 3.5, when we were restricted to intervals of ν-measure 1/2, and found a maximal invariant set W r νn ([0.5001, 0.8489]) = 0.6821, which is much lower than any set that we find by Methods 1-4.

Examples
We demonstrate the ideas from the previous sections using a more complex numerical example.For clarity, we first outline the steps to carry out this process in the following algorithm.The code to execute Steps 3 onwards is available on the authors' websites, and Steps 1 and 2 are easily carried out using MATLAB and the purpose-built software GAIO [12].Algorithm 1.
1. Partition the state space X r into connected sets {B 1 , B 2 , . . ., B n }.This can be achieved relatively easily using eg.GAIO via MATLAB.GAIO [12] integrates the partitioning process with the construction of the transition matrix.
2. Construct the transition matrix corresponding to the closed system, P r = (P r . The matrix P r can be quickly estimated using images of sample points within each box, see eg. [12].
3. Compute the vector p r as the leading left eigenvector of P r , normalised to sum to 1.This is fast using eg.MATLAB as P r is sparse.
5. Form Pij = p r j Pji /p r i and R = ( P + P )/2, and compute the second largest eigenvalue λ R 2 and corresponding right eigenvector v of R.

Return max{W r
νn (A 1 (c 1 )), W r νn (A 2 (c 2 ))} as the maximal almost-invariance ratio found for the open system and set A opt to be whichever of A 1 (c 1 ) or A 2 (c 2 ) achieves this maximum.
We will now apply Algorithm 1 to a periodically forced two-dimensional system of ODEs defined on X = [0, 2] × [0, 1] by: ẋ = −πA sin(πf (x, t)) cos(πy) (5.1a) Lebesgue measure is preserved by the corresponding flow.This flow, known as the "Double Gyre" has been studied in several recent papers including [42,29].Following [29], we set A = 0.25, ǫ = 0.25 and ω = 2π so that the period of the flow equals 1.We divide the state space X into a 2 8 × 2 7 grid, so that n = 2 15 .The matrix P is then formed by considering ).After removing these holes we have n = 17198 boxes covering the restricted state space X r .The leading eigenvalues of P r are λ r 1 = 0.9740 and λ r 2 = 0.8679.The leading left eigenvector p r is then computed (see Figure 9), and we use this to remove all boxes with ν n -measure of less than 10 −9 .This leaves us n = 12843 boxes.We recompute P r , and construct P , R and v 2 as described in Steps 4-5 (see Figure 10).We then perform the thresholding process described in Steps 6-7 with the second right eigenvector of R. The invariance ratios W r νn (A 1 (c)) and W r νn (A 2 (c)) obtained for different threshold values c are plotted in Figure 11.The maximal invariance ratio for a set of ν n -measure less than 1/2 is found when c * = 0.001805, and then we have This set has ν n -measure 0.4565 and an invariance ratio W r νn (A 1 ) = 0.9305, and is shown in Figure 12.For comparison, we also repeat Steps 6-7 using the second left eigenvector of P r , and find a set A opt , shown in Figure 13.The results are compared in Table 2; we once again find a higher invariance ratio by thresholding v 2 than any other candidate eigenvector.
We note that the separation seen in Figure 12 is similar to the separation seen in the closed system (Figure 16 in [29]), despite the significantly different asymptotic behaviour; the conditionally invariant measure shown in Figure 9 is significantly different to the areapreservation of the closed system.As the separation seen in [29] appears to be strongly related to the stable and unstable manifolds of the hyperbolic periodic points on the top and bottom of the rectangle, and at least part of these manifolds are contained in the support of the conditionally invariant measure in Figure 9, it may be that transport related to lobe dynamics still plays a dominant transmission role in the open system; we refer the reader to [29] for further discussion of the lobes.

Conclusion
The determination of almost-invariant or metastable regions in dynamical systems provides important information on the transport of mass at a macroscopic level.The most successful methods to study almost-invariance in closed systems are based on the Perron-Frobenius operator P. While considerable research has been done on almost-invariance in closed systems, little, if any attention has been paid to open systems, where trajectories terminate when entering a hole, and the time-asymptotic dynamics occurs on a survivor set.We have introduced the notion of almost-invariance in open systems and developed a new closing operation to construct a Markov operator P from the conditional Perron-Frobenius operator P r corresponding to the open system.Existing techniques for almost-invariant sets, such as those in [27], may now be applied to this Markov operator.We showed that the invariance ratios of sets according to the (i) the conditional Perron-Frobenius operator P r and (ii) the Markov operator P , differ by an amount equal to the conditionally invariant measure of the set in question intersected with the hole.We further developed rigorous upper and lower bounds for the maximum possible invariance ratios of discretised sets based on the second eigenvalue of a matrix constructed from the closing operation.
The new methodology was demonstrated to be superior to naive applications of existing approaches in one-dimensional and two-dimensional examples.One drawback with the approach presented here is the non-sparsity of the matrix resulting from the closing operation.Although one can still compute with tens of thousands of boxes (eg.12843 boxes in the as required.Secondly, we show that f r is a fixed point of P. We note which shows that f r is a fixed point of P. Lastly, it remains to relate the spectra of P r and P. For all g ∈ L 1 (ℓ r ) such that P r g = λ r g with |λ r | < 1, we have X r g dℓ r = 0, and therefore Pg = P r g.Hence, Pg = λ r g.Proof of Lemma 4.4.As Q = π r P = π r (P r + τ • f r ), our main task is to construct the matrix representation for π r (τ (g) • f r ) when g ∈ span{½ B i : i ∈ I r }.We write g = i∈I r g i ½ B i ,

Figure 6 :
Figure 6: The invariance ratios for the closed system and the open system, for all intervals of invariant measure 1/2 for the 1/16-shifted doubling map.The dashed lines in the second panel indicates the preimage of the hole.The graph of T is overlaid onto the second panel.

Figure 7 :
Figure7: The difference Wν (A) − W r νn (A) for intervals of ν-measure 1/2 for the 1/16-shifted doubling map.The difference is zero when the interval does not include the preimage of the hole (which is demarcated with dashed lines), and is bounded from above by 1−λ r n,1 = 0.0361.

Method 2
Threshold the second right eigenvector of R (Figures8c, 8d

) 4 . 3 2 ⊥
We use the right eigenvectors of R and R because one has v right p (recall p is the leading left eigenvector) by a standard result of linear algebra.Thus we think of the entries in v right 2 as relaxations of ±1, where the value +1 corresponds to membership in one almost-invariant set and −1 corresponds to membership in the complementary almost-invariant set.

Figure 8 :
Figure 8: The graphs in the left column depict the eigenvector that was used for thresholding, together with the optimal set A opt The graphs in the right hand column show the threholding process.The red lines corresponds to W r νn (A 1 (c)) and the blue dashed line to W r νn (A 2 (c)).The black horizontal line represents the value for which ν n (A 1 (c)) = ν n (A 2 (c)) = 1/2.The value of c * corresponding to the W νn (A opt ) is indicated.

Figure 9 :
Figure 9: The conditionally invariant measure of the double gyre defined on X r = X \ (H 1 ∪ H 2 ).

Figure 13 :
Figure13: For the periodically forced double gyre defined on X r , the partition into optimally almost-invariant sets found by thresholding with the second left eigenvector of P r .The set with the highest invariance ratio of 0.8527 is shown in orange, and the complementary set is shown in dark brown.

Table 1 :
Summary of results for α-shift map.Method Vector used for thresholding W νn (A opt ) ν n (A opt )