Support vector machines for learning reactive islands

We develop a machine learning framework that can be applied to data sets derived from the trajectories of Hamilton’s equations. The goal is to learn the phase space structures that play the governing role for phase space transport relevant to particular applications. Our focus is on learning reactive islands in two degrees-of-freedom Hamiltonian systems. Reactive islands are constructed from the stable and unstable manifolds of unstable periodic orbits and play the role of quantifying transition dynamics. We show that support vector machines (SVM) is an appropriate machine learning framework for this purpose as it provides an approach for ﬁnding the boundaries between qualitatively distinct dynamical behaviors, which is in the spirit of the phase space transport framework. We show how our method allows us to ﬁnd reactive islands directly in the sense that we do not have to ﬁrst compute unstable periodic orbits and their stable and unstable manifolds. We apply our approach to the H´enon-Heiles Hamiltonian system, which is a benchmark system in the dynamical systems community. We discuss diﬀerent sampling and learning approaches and their advantages and disadvantages.

Prediction and control of transition between potential wells across index-one saddles are of importance in a diverse array of non-linear systems in physics, chemistry, and engineering.Dynamical systems theory provides the framework for understanding transition dynamics in these systems.The relevant phase space structures are the stable and unstable invariant manifolds of the hyperbolic periodic orbit associated with the index-one saddle.Thus, to determine the fate of an initial condition in the well without solving the non-linear equations, one has to check whether the initial condition lies in a region bounded by the globalized stable (escape from the well) or unstable (entering the well) invariant manifold.The region that bounds the initial condition is given by the intersection of the invariant manifold with a two-dimensional section and called the reactive island, borrowing terminology from chemical reaction dynamics.In this paper, we develop and verify a trajectory-based framework using support vector machines, by applying it to the Hénon-Heiles Hamiltonian, for learning the reactive islands.Our results show that support vector machines are an ideal data-driven framework for learning the geometry of phase space structures.The approaches developed here are robust to changes in system parameters and geometry of the reactive islands.

I. INTRODUCTION
The goal of this paper is to develop a machine learning framework that identifies the phase space structures governing phase space transport in data sets constructed from trajectories of Hamilton's equations.Our focus in this paper will be on two degrees-of-freedom Hamiltonian systems.
The Hamiltonian function typically has the form of the sum of a kinetic energy, a function of the momentum variables and the potential energy, a function of the position variables.For Hamilton's equations in canonical form each momentum variable is canonically conjugate to one position variable.Hence there are an equal number of momentum and position variables.
The space of both the momentum and position variables is referred to as the phase space and the space of only the position variables is referred to as the configuration space.The number of degrees of freedom of a Hamiltonian system is the number of configuration space variables.
Hamilton's equations describe dynamics in phase space.Nevertheless, the topography of the potential energy surface plays an important role in defining the reaction dynamics problem for a specific Hamiltonian system.In particular, potential wells are identified as reactants and products.Wells are separated by saddle points and the configuration space picture of reaction dynamics is that of a trajectory evolving between wells by "crossing" the saddle point.
Index-one saddle points on a potential energy surface (PES) are the "seeds" for the phase space structures from which the theory of reactive islands is constructed.Conley 1,2 was the first to analyze the phase space geometry and dynamics near an index-one saddle in two degrees-of-freedom (DoF) Hamiltonian systems in his studies of the circular restricted three body problem.
The Lyapunov subcenter theorem 3,4 is fundamental for passing from an equilibrium point (the index-one saddle) to dynamical behavior.It states that for a range of energies above the energy of the index-one saddle (the exact range is not given in the theorem) there exists an hyperbolic periodic orbit having two-dimensional stable and unstable manifolds.
In the three-dimensional energy surface, stable and unstable manifolds have the geometry of cylinders ("tubes") and these stable and unstable cylinders mediate transport across the region of the index-one saddle.
The cylindrical structure of the stable and unstable manifolds, together with their invariance, implies that trajectories starting within the tubes must remain in the tubes for all positive and negative time.All trajectories inside the stable cylinder approach the hyperbolic periodic orbit in positive time, pass through the region bounded by the periodic orbit, and then exit the region through the unstable cylinder.The stable cylinder (and the unstable cylinder) has two "branches" joined at the periodic orbit, one emanating to each "side" of the orbit.
In the context of escape from a potential well, escaping trajectories must lie in a branch of the stable cylinder.Trajectories starting in the branch of the stable cylinder lying in the potential well correspond to forward escape trajectories, i.e. trajectories that escape in forward time.Trajectories starting in the branch of the stable cylinder lying outside the potential well correspond to backward escape trajectories, i.e. trajectories that are captured in the well in forward time.After escape from/capture in the well, trajectories are guided by an unstable cylinder.This is how stable and unstable cylinders govern the dynamics in phase space.
We describe a lower dimensional technique for probing the geometry of the cylinders.In the potential well we construct a two dimensional Poincaré section, i.e. a two dimensional surface where the Hamiltonian vector field is everywhere transverse to the surface.The Poincaré section is constructed in such a way that the stable (and unstable) cylinder intersects it in a topological circle.The region bounded by this topological circle is referred to as a "reactive island", where 'reactive' refers to the occurence of a chemical reaction, an analogue of the escape from a potential well considered here.The Poincaré map of the Poincaré section into itself is the map that associates to a point its first return to the Poincaré section under the flow generated by the Hamiltonian vector field.The inverse of this Poincaré map of the Poincaré section into itself is the map that associates to a point its first return to the Poincaré section under negative time, i.e. the point "where it came from".
We consider the preimages of this reactive island by considering its evolution backwards in time under the inverse of the Poincaré map.In this way one obtains a reactive island on the Poincaré section that returns to the original reactive island, and then escapes for positive time evolution.By repeating this construction we obtain a sequence of reactive islands, ordered in time, which sequentially map to each other before reacting.Possible pathological intersections of the cylinders with the Poincaré section can occur, and are discussed in 5 .A theory of reaction dynamics for two DoF systems based on the geometry of these stable and unstable cylinders was developed in the late 1980's and 90's in [5][6][7][8][9][10][11][12] , although some of the ideas appeared in earlier work 13,14 , and it goes by the name of "reactive island theory".
Our goal is to consider a data set consisting of points on the energy surface that are labeled based on the evolution of the corresponding trajectory and "learn" which points lead to a escape from a potential well and which do not.We have chosen to use support vector machines (SVM), a class of supervised learning classification and regression algorithms [15][16][17] , which has been introduced to nonlinear dynamical systems by Ref. 18.We employ the classification algorithms, which define a decision function by determining the boundary between different classes of data.To be precise, SVM identifies a subset of the data set referred to as "support vectors" to calculate the boundary for which the distance to data in both classes is locally maximal.In our case, the classes in the classification correspond to initial conditions leading to "qualitatively different" dynamical behaviour.The learned boundary between escaping and non-escaping trajectories in phase space consists of the invariant manifolds discussed above and thereby SVM enables us to determine the geometrical structures in phase space governing reaction dynamics.
The reasons for choosing support vector machine (SVM) to identify the phase space structures are: (i) Nonlinear kernels in SVM provide means to approximate curves such as reactive islands which form nonlinear boundaries between reactive and non-reactive trajectories.We would like to note that methods for nonlinear clustering 19 and other kernel methods 20 are also candidates for developing similar approaches, but beyond the scope of this study.(ii) We design our approach with extensions to higher-dimensional applications in mind.SVM is known to work well even with small amounts of data, therefore our approach is computationally better suited than existing methods for high-dimensional systems and systems where integrating trajectories is expensive.This makes SVM suitable for the reactive island theory of three DoF 21 and system-bath models of isomerization 22 which has been developed recently and supports the case for a computationally efficient approach to phase space structures presented here.
The Hamiltonian system that we will use to benchmark our approach will be the Hénon-Heiles system 23 .This is a two degree of freedom Hamiltonian system that serves as a paradigm for understanding complex dynamics in a variety of settings.As a function of the energy, it can display dynamical behavior that numerically appears to be "near integrable" to completely chaotic.Moreover, this system has three index one saddles that define three distinct reaction channels having the geometric structure discussed above.The geometry of reaction dynamics geometry for a similar system has been analyzed and discussed in Ref. 24 .This paper is outlined as follows.In Section II we describe the two degrees-of-freedom Hénon-Heiles Hamiltonian system and the nature of reaction dynamics in the context of this model.In Section III we describe the machine learning technique of support vector machines (SVM) and the approach of active learning.We describe how trajectories are constitute into data sets and the use of Lagrangian descriptors as a trajectory diagnostic.In Section IV we describe our results, and in Section V we present our conclusions and the outlook for further research.

II. H ÉNON-HEILES SYSTEM AND REACTION DYNAMICS
We use the Hénon-Heiles Hamiltonian with three index-one saddles in bottlenecks through which trajectories can escape.This escape can be interpreted as a reaction by crossing the potential energy barrier.This model system and its high dimensional analog has been studied in great detail in nonlinear dynamical systems, statistical mechanics, for developing molecular simulation algorithms [25][26][27][28][29] .In this study, we will define the three exits as follows: the entering via the top index-one saddle corresponds to the formation of a molecule complex by combination of two atoms or molecules, while the left and right exits correspond to dissociation of the molecule into two different products with structural symmetry.
The Hamiltonian is given by where all the parameters are set to 1.0 to be comparable with the known results in the literature 23,30,31 .The vector field is given by and the equilibrium points are (0, 0, 0, 0) , 0, The eigenvalues of the linearized system in Appendix A1 at the equilibrium points gives the origin is a center × center equilibrium and the three remaining are saddle × center equilibrium points.
We illustrate the dynamics given by the vector field (2) by showing four sample trajectories in the Fig. 1 at same energy starting at the center × center equilibrium point with different p x momenta and p y > 0. In Fig. 1 (a,c), p x coordinates only differ in the second significant digit, and yet the escape from the potential well occurs via different bottlenecks and at different times.However, in Fig. 1 (b,d), a similar difference in p x coordinates only changes the escape time.This illustrates the challenge of sampling an ensemble of escape (transition) trajectories.Furthermore, the challenge in identifying different timescales of escape trajectories is due to the large variation in time to escape from the potential well when there is merely a small difference in p x .We will next show the underlying phase space structure of such escape behavior and then show the use of trajectory data in identifying the structure.Finally, it is important to note that for the total energies that we consider the energy surface is unbounded.This implies that notions of recurrence from standard ergodic theory, such as Poincaré recurrence and ergodicity, do not apply because escaping trajectories become unbounded and never return.
The trajectory behavior shown in Fig. 1 is mediated by the stable manifolds (cylindrical geometry or tubes) of the hyperbolic periodic orbits associated with the index-one saddles.
For the two degrees-of-freedom Hamiltonian system, tube manifolds can be computed and visualized in their complete geometry as shown in Fig. 2 (a) for the index-one saddle at 0, ω 2 y δ , 0, 0 .The stable manifolds associated with the three saddles are projected on the configuration space (x, y) are shown in Fig. 2 (b) (In Appendix.A, Fig. 9 shows all the tube manifolds in 3D).In this article, we compute the stable manifolds using the procedure described in the Appendix A to obtain the segments which intersect a section with y = 0 and p y > 0 shown in Fig. 2(b).Given such a section of the three dimensional energy surface, the first order reactive islands of escape are defined as the last intersection of the stable manifolds with the section.By last intersection, we mean the trajectories in forward time intersect the section and then leave the potential well without returning to the section and also referred to as the imminent escape.

III. SUPPORT VECTOR MACHINES AND ACTIVE LEARNING
The classification algorithms for SVM [15][16][17]  complementary area that does not lead to imminent escape.The exact boundaries between these areas are the reactive islands formed by stable and unstable invariant manifolds of hyperbolic periodic orbits discussed above.A SVM classification algorithm, also referred to as support vector classifier (SVC), therefore approximates reactive islands in this setting.
Similarly to 32, we use the scikit-learn 33 implementation of SVM 34 .The implementation can be used with various kernels, of which the radial basis function kernel is best suited to approximate reactive islands in the Hénon-Heiles system, which are topological circles.
With this kernel, a previously unseen data point P is predicted to belong to a class using the decision function where γ > 0 controls the width of the Gaussian, l i = ±1 are class labels of training data P i and C ≥ α i ≥ 0 are weights calculated by the algorithm, of which only a number of the weights α i is non-zero.These weights correspond to the training data P i called support vectors.The weights α i are calculated by SVC such that the distance between the predicted boundary and the closest points P i of every class is maximised, as illustrated in Fig. 3.While it is possible to apply SVM to a fixed training data set, such as a regular grid, the accuracy of the resulting decision boundary will be limited by the amount of training data and its spacing.A significantly less data-intensive approach is offered by active learning 35,36 , where the 'learner' biases its sampling based on information obtained from previous samples.
To do this, we start SVM with a coarse grid of data points and iteratively add data points in the proximity of the decision boundary and re-run SVM.This allows the algorithm to explore the intricate structures usually formed by invariant manifolds in systems describing chemical reaction dynamics.At every iteration we randomly add one data point near 10 randomly selected support vectors.The point is added using the multivariate normal distribution N (P, I), where P is the support vector and I the identity matrix.
We would like to point out the importance of the precise problem formulation.Homoclinic and heteroclinic intersections of invariant manifolds lead to fractal structures, that is a fractal boundary between classes of dynamics.There is no known way to resolve fractal structures with finite precision and a finite number of data points.Thus approximating the boundary using fixed-width Gaussians is bound to fail.In many systems it is possible to avoid these fractal structures by carefully selecting a surface of initial conditions and studying the dynamics under the corresponding return (Poincaré) map.

IV. RESULTS AND DISCUSSION
In the figures below, the cyan dots denote the support vectors used by the classifier to learn the decision boundary and the black dashed line denotes the learned decision boundary.
The true reactive islands are shown as red, green, and blue curves.
Fixed training data.We sample initial conditions on the two-dimensional section of the three dimensional energy surface, H(x, y, p x , p y ) = E, with y = y c and p y > 0. In this study, we present the results for two sections at y c = 0, −0.25 to denote distinct sections at the location of the well and near the bottleneck.We also show the results for E = 0.17 − 0.20 in increments of 0.1 to illustrate the approach for various imbalance (ratio of reactive to non-reactive trajectories) in the training data corresponding to different excess energies.We generate a grid of initial conditions 100 × 100 and sample the y−momenta using the fixed energy condition.Then, we run trajectories with the initial conditions for a prediction time horizon of t = 30 time units.Then, we classify the escape trajectories as reactive through bottleneck 1 or 2 or 3 if they cross the line x = −1.25, or x = 1.25, or y = 1.25, respectively.
If an initial condition does not satisfy any of the above conditions for the chosen time interval, we label it as non-reactive denoted by 0. Thus, we obtain a multi-label training dataset with two feature vectors, x, p x coordinates to discover the distribution of reactive trajectories on the two-dimensional section corresponding to the two coordinates.We use a smaller interval for the Gaussian radial basis function parameters, C ∈ {1e2, 1e3, 1e4, 1e5} and γ ∈ {10, 1e2, 1e3, 1e4}, to perform a grid search for high accuracy along with 5-fold cross-validation during the training.
The decision boundaries as learned by the SVC trained using fixed size training dataset and the reactive islands obtained from the direct computation of the tube manifolds are compared for verification in Fig. 5.The learned decision boundaries track the true reactive islands to a high accuracy and the classification accuracy is above 99% for all the energies considered here.In fact, for the low energy case (E = 0.17, which is 0.003 above the energy of the saddle, the accuracy is similar to the highest energy case, 0.033.These two energies represent two extreme training data sets as the fraction of reactive trajectories increases with total energy.Thus, for low energy case a uniform sampling is bound to give small number of reactive trajectories and vice versa for high energy case.Even though this can be corrected using a non-uniform sampling for low energy, we show that learning the reactive islands leads to high accuracy because it is the fundamental phase space structure underlying the reactive trajectories.
Active learning.In this approach, we developed an iterative method for active learning The decision boundaries as learned by the SVC using an active learning approach and the reactive islands obtained from the direct computation of the tube manifolds are compared for verification in Fig. 6.
Trajectory geometry enabled learning.In this approach, we use a positive scalar quantity for encoding the geometry of a trajectory called the Lagrangian descriptor (LD) 37,38 to generate a new feature for training a SVC.Lagrangian descriptors have been shown to detect phase space structures that mediate transition dynamics in general non-linear dynamical systems (see Ref. 39 and references therein).Furthermore, as a trajectory based diagnostic of the phase space, it can be computed on-the-fly with the trajectory.This gives it a two-fold merit as a feature: (i) encodes the geometry of the trajectory, thus incorporates the phase space perspective and (ii) efficient computation along with trajectory generation.
We briefly describe the method of Lagrangian descriptors which reveals regions with qualitatively distinct dynamical behavior by showing the intersection of the invariant manifolds with the two dimensional section.For a general time-dependent dynamical system given by where the vector field f (x, t) is assumed to be sufficiently smooth both in space and time.The vector field f can be prescribed by an analytical model or given from numerical simulations as a discrete spatio-temporal data set.For instance, the vector field could represent the velocity field of oceanic or atmospheric currents obtained from satellite measurements or from the numerical solution of geophysical models.For any initial condition x(t 0 ) = x 0 , the system of first order nonlinear differential equations given in Eqn.
(5) has a unique solution represented by the trajectory that starts from that initial point x 0 at time t 0 .
In this study, we adopt the LD definition where f k is the k−the component of the vector field, Eqn. ( 5) and use p = 1/2.We note that the integral can be split into its forward and backward time parts to detect the intersection of stable and unstable manifolds separately.This relates to finding the escape and entry channels into the potential well.In this study, we keep the forward part of the integral given by Although this definition of LD does not have an intuitive physical interpretation as that of the arclength definition 38 , it allows for a rigorous proof that the "singular features" (nondifferentiable points) in the LD contour map identify intersections with stable and unstable invariant manifolds 40 .Another important aspect of what is known in LD literature as the p-(quasi)norm is that degrees of freedom with relevance in escape/transition (reaction) dynamics can be decomposed and computed.This definition was used to show that the method can be used to successfully detect NHIMs and their stable and unstable manifolds in Hénon-Heiles Hamiltonian 24,31 .For this system, where both fixed (or variable) integration time is used, it has also been shown that the LD scalar field attains a minimum (or maximum) value along with singularity at the intersections of the stable and unstable manifolds, and given by where W s (x 0 , t 0 ) are the stable manifolds calculated at time t 0 and argmin denotes the phase space coordinates on the two dimensional section that minimize the scalar field, L f p (x 0 , t 0 , τ ), over the integration time, τ .Thus, the scalar field plotted as a contour map identifies the intersection of the stable manifold with a two dimensional section.This ability of LD contour map to partition trajectories with different phase space geometry is shown in Fig. 7(a-c) as values of LD inside the reactive islands are close to constant.
We construct the training data with three features -x, p x , M 0.5 (τ ) -and a fixed size dataset.Then, we implement a SVC as in the Fixed training data approach with the parameters for the Gaussian radial basis function C ∈ {1e2, 1e3, 1e4, 1e5} and γ ∈ {10, 1e2, 1e3, 1e4}.
In Fig. 8 (d-i), we show the reactive islands along with the predictions of a SVC trained using the trajectory geometry given by LD as a feature and with a training data size of 10000 points.We note that when such a data set is used the support vectors used by the model increases around the boundary.However, this approach encodes the geometry of a trajectory in phase space and hence leads to robust classification as the total energy parameter is increased.

V. CONCLUSIONS AND OUTLOOK
On a transverse two dimensional section of the energy surface, reactive islands are the intersections of stable and unstable invariant manifolds of hyperbolic periodic orbits.Thus, the one dimensional boundaries of the reactive islands separate transition and non-transition trajectories.In this article, we presented three support vector classifier approaches for learning the reactive islands: fixed dataset training, active learning, and trajectory geometry enabled training.The advantages of our approach are as follows: (a) avoiding the need to compute hyperbolic periodic orbits and the associated invariant manifolds, in favour of finding the reactive islands directly on a surface of section as the boundary between classes of qualitatively different dynamics, (b) minimising computational cost of trajectory calculations by sampling the section near a boundary learned from a coarser sampling (c) using trajectory geometry as a dynamical (phase space) feature in the training data and compressing the high dimensional trajectory into a smaller feature set.Inheriting low data requirements from SVM, our approach is designed to work well for systems where integrating trajectories is expensive and is expected to generalise well systems with more than two degrees of freedom.
Our work intends to simplify the process of finding reactive islands, making it easier to generalise and more accessible to a wider scientific audience.Future work in this direction will involve the application to a model of chemical reaction and examples of high dimensional phase space of a system-bath model.

VII. DATA AVAILABILITY
The code that support the findings of this study are openly available on the corresponding author's GitHub repository.The data used in training the machine learning model is available from the corresponding author upon reasonable request.
Appendix A: Hyperbolic periodic orbit and invariant manifolds

The Linearized Hamiltonian System
The linearized equation of motion around the index-one saddle equilibria with coordinates (x e , y e , 0, 0) is given by expanding the terms of the Hamiltonian (1) and keeping the quadratic terms.After making a coordinate to (x e , y e , 0, 0) as the origin, the quadratic Hamiltonian function gives the linear system at the equilibrium point The linearized dynamics near the index-one saddle equilibrium points extend to the full nonlinear system due to Moser's generalization of Lyapunov's theorem.
Thus, the eigenvectors corresponding to , β = ±λ, ±iω, can be written as respectively.Thus, the general solution of the linear system near the saddle equilibrium point is given by with A 1 , A 2 being real and B = B 1 + iB 2 being complex.

Computing the hyperbolic periodic orbit and associated tube manifolds at the index-one saddle
For discussing the geometry, we call the equilibrium with positive y-coordinate x eq,top and negative y-coordinate x eq,left and x eq,right .
Select appropriate energy above the critical value -For computation of manifolds that act as boundary between the escape and non-escape trajectories, we select the total energy, E, above the energy of the index-one saddle, E s , and thus the excess energy Differential correction and numerical continuation -We present a procedure which computes the hyperbolic periodic orbits associated with an index-one saddle using order correction δy 0 as which is iterated until |p y 1 | = |δp y 1 | < for some tolerance , since we want the final periodic orbit to be of the form This procedure yields an accurate initial condition for a periodic orbit of small amplitude A x ≈ 10 −4 , since our initial guess is based on the linearization at the equilibrium point.
Numerical continuation to periodic orbit at arbitrary energy.-Theprocedure described above yields an accurate initial condition for a periodic orbit from a single initial guess.
If our initial guess came from the linear approximation near the equilibrium point, from Eqn. (A7), it has been observed numerically that we can only use this procedure for small amplitude, of order 10 −4 , periodic orbits around x eq,bot .This small amplitude correspond to small excess energy, typically of the order 10 −2 , and if we want to compute the periodic orbit of arbitrarily large amplitude, we resort to numerical continuation for generating a family which reaches the appropriate total energy.This is done using two nearby periodic orbits of small amplitude to obtain initial guess for the next periodic orbit and performing differential correction to this guess.To this end, we proceed as follows.Suppose we find two small nearby periodic orbit initial conditions, x(1) 0 and x(2) 0 , correct to within the tolerance d tol , using the differential correction procedure described above.We can generate a family of periodic orbits with successively increasing amplitudes around xeq,bot in the following way.

Let ∆ =
x(2) 0 − x(1) 0 = [∆x 0 , ∆y 0 , 0, 0] T (A19) A linear extrapolation to an initial guess of slightly larger amplitude, x(3) 0 is given by x(3) 0,g = x(2) 0 + ∆x 0 ), (y 0 + ∆y 0 ), 0, 0 0 , 0, 0 T (A20) so small that the time of flight becomes too large due to asymptotic nature of the stable and unstable manifolds.Ref. 41 suggests typical values of > 0 corresponding to nondimensional position displacements of magnitude around 10 −6 .By numerically integrating the unstable vector forwards in time, using both and − , for the forward and backward branches respectively, we generate trajectories shadowing the two branches, W u + and W u − , of the unstable manifold of the periodic orbit.Similarly, by integrating the stable vector backwards in time, using both and − , for forward and backward branch respectively, we generate trajectories shadowing the stable manifold, W s +,− .For the manifold at X(t), one can simply use the state transition matrix to transport the eigenvectors from X 0 to X(t) as X s (X(t)) = Φ(t, 0)X s (X 0 ) (A24) It is to be noted that since the state transition matrix does not preserve the norm, the resulting vector must be normalized.The globalized invariant manifolds associated with index-one saddles are known as Conley-McGehee tubes 46 .These tubes form the skeleton of transition dynamics by acting as conduits for the states inside them to travel between potential wells.
The computation of codimension-1 separatrix associated with the hyperbolic periodic orbit around a index-one saddle begins with the linearized equations of motion.This is obtained after a coordinate transformation to the saddle equilibrium point and Taylor expansion of the equations of motion.Keeping the first order terms in this expansion, we obtain the eigenvalues and eigenvectors of the linearized system.The eigenvectors corresponding to the center direction provide the starting guess for computing the hyperbolic periodic orbits for small excess energy, ∆E << 1, above the saddle's energy.This iterative procedure performs small correction to the starting guess based on the terminal condition of the periodic orbit until a desired tolerance is satisfied.This procedure is known as differential correction and generates hyperbolic periodic orbits for small excess energy.Next, a numerical continuation is implemented to follow the small energy (amplitude) periodic orbits out to high excess energies.

FIG. 1 .
FIG.1.Four sample trajectories initialized on the section y = 0, p y > 0 with p x = 0.516, 0.07, 0.526, 0.08 in (a-d), respectively, projected on the configuration space.The total energy E = 0.17 is slightly above the energy of the index one saddles and escape times T E are shown on each plot.
FIG. 2. (a) Cylindrical (or tube) manifolds, stable in green and unstable in red, of the hyperbolic periodic orbits associated with the top index-one saddle in the Hénon-Heiles Hamiltonian.The energy of the hyperbolic periodic orbit and the invariant manifolds are at the total energy, E = 0.17 and mediate the trajectories that escape via top saddle as shown in Fig. 1(b,d).(b) Stable manifolds projected on the configuration space reveal the geometry of imminent escape from the potential well via the three bottlenecks.Only the segment of the stable manifolds from the hyperbolic periodic orbits to the intersection with the Poincaré section (shown as a black line) is shown for the energy, E = 0.19.

FIG. 3 .
FIG. 3.An illustration of a decision boundary (black) between two classes of data (blue and orange)calculated by SVC with radial basis function kernel.The distance between the boundary and the closest points P i of every class, in this case the support vectors, is highlighted in green.

FIG. 4 .
FIG. 4. Colormap showing accuracy for different combination of radial basis function parameters, (C, γ).The accuracy is obtained using a 5-fold cross validation over the grid of values for C and γ.The pair of value which gives maximum accuracy is chosen for training the support vector classifier.

25 FIG. 5 .
FIG. 5. Fixed training data set.Reactive islands identified by the support vector classifier trained using fixed size data set shown as dashed curves.The overlayed continuous curve is obtained using direct computation of tube manifolds at energy E = 0.17, 0.18, 0.19, 0.20 in first, second, third, fourth column, respectively.The magenta curve denotes the intersection of the energy surface with the two dimensional section.The cyan dots denote the support vectors used by the classifier in learning the reactive islands as decision surfaces.Two sections with (x, p x ) coordinates are shown in top and bottom rows: (a-d) y c = 0 (e-h) y c = −0.25 with p y > 0.

25 FIG. 6 .
FIG. 6. Active learning Reactive islands identified by the support vector classifier trained using data generated near the coarse boundaries shown as dashed curves.The overlayed continuous curve is obtained using direct computation of tube manifolds at energy E = 0.17, 0.18, 0.19, 0.20 in first, second, third, fourth column, respectively.The magenta curve denotes the intersection of the energy surface with the two dimensional section.The cyan dots denote the support vectors used by the classifier in learning the reactive islands as decision surfaces.Two sections with (x, p x ) coordinates are shown in top and bottom rows: (a-d) y c = 0 (e-h) y c = −0.25 with p y > 0.

FIG. 7 .
FIG. 7. Reactive islands at E = 0.18, 0.20 are shown as red, green, and blue curves, respectively, and the LD values corresponding to the initial conditions sampled on y c = 0 with p y > 0 section.These initial conditions along with the LD value computed for τ = 30 or until a trajectory escapes is used as the training data.

25 FIG. 8 .
FIG. 8. Trajectory geometry enabled learning Shows the reactive islands, support vectors (as dots), and predictions (as cross) of the trajectory geometry trained support vector classifier.