Koopman Operator and Phase Space Partition of Chaotic Maps
aa r X i v : . [ n li n . C D ] J u l Koopman Operator and Phase Space Partition of Chaotic Maps
Cong Zhang and Yueheng Lan
1, 2, ∗ School of Science, Beijing University of Posts and Telecommunications, Beijing 100876, China State Key Lab of Information Photonics and Optical Communications,Beijing University of Posts and Telecommunications, Beijing 100876, China (Dated: July 23, 2020)Koopman operator describes evolution of observables in the phase space, which could be used to extractcharacteristic dynamical features of a nonlinear system. Here, we show that it is possible to carry out interestingsymbolic partitions based on properly constructed eigenfunctions of the operator for chaotic maps. The partitionboundaries are the extrema of these eigenfunctions, the accuracy of which is improved by including more basisfunctions in the numerical computation. The validity of this scheme is demonstrated in well-known 1-d and 2-dmaps. It seems no obstacle to extend the computation to nonlinear systems of high dimensions, which providesa possible way of dissecting complex dynamics.
PACS numbers: 05.45.-a, 02.30.Tb, 05.45.Ac, 02.70.WzKeywords: Koopman Operator, Phase Space Partition, Symbolic Dynamics, Dynamical Model, Chaotic Map
I. INTRODUCTION
For complex systems, due to the lack of understanding ofunderlying physical principles or generic complexity, it is im-possible to make an effective description or even pin down themajor characteristics of its dynamics. With the rapid develop-ment of data collecting and storing technology, however, it isindeed feasible to obtain large amount of data through large-scale simulation, detailed experiments or high-resolution ob-servation. Although these data may be handled with modernintelligent "black box", e.g. , various neural networks [1], adeeper understanding of key dynamical features of the sys-tem remains desperately desirable and thus of essential im-portance.Dynamical systems theory often studies orbit structures inphase space and their changes under parameter variations. Inlinear or integrable systems, all the orbits are arranged regu-larly and qualitative features of the system may be estimatedwith very high accuracy. However, when the dynamics be-comes chaotic, which is the case for most nonlinear systemsin high dimensions, traditional analytic approaches often donot apply [2]. Even numerical computation becomes unreli-able [3] in these systems, where the orbit structure becomesextremely complicated and sensitive to initial condition or ex-ternal perturbation. It seems necessay to introduce a statisticalassumption in this context like what has been done in statisti-cal mechanics [4]. However, the lack of a variational princi-ple and non-stationarity in general requires that it should startfrom first principles. Otherwise important dynamical featuresmay be missed in the treatment of nonlinear systems far fromequilibrium. This concern has been raised long time ago in aneffort to build rigorous framework of statistical physics, whereLiouville equation is viewed as a starting point by many [5].In this type of approach, the focus is shifted from individ-ual orbits to the evolution of specific functions defined in thephase space. In Hamiltonian systems, this job is done by ∗ [email protected] the well-known Liouville operator while in general dynami-cal systems, Koopman operator is the analogue.The Koopman operator was put forward by B. O. Koopmanin 1931 [6], which describes evolution of functions defined inthe phase space of a dynamical system and is closely relatedto Liouville operator [7, 8]. The eigenvalues and eigenfunc-tions of this linear operator capture global features of systemdynamics [9], which were first used by Mezi´c et al for ergodicpartitions in systems with heterogeneous dynamics [10], andfor the possibility of simplifying high-dimensional dynamicalsystems [11]. In recent years, the dynamical mode decom-position (DMD) algorithm based on Koopman operator andits improved versions have been extensively tested on differ-ent occasions including the power system [12, 13], buildingenergy efficiency [14, 15] and fluid systems [16–18]. Therobustness of the DMD algorithm is discussed in [19] andthe fine structure of the spectrum could in principle be com-puted with rigorous convergence guarantee from measureddata [20]. In practice, Hankel type of matrices is easy to con-struct for a convenient computation of the spectra by properlyarranging the data series [21]. Further exploration reveals aconnection of the linearizing transformation to the spectrumof Koopman operator [22], which provides a link to anotherpowerful tool in the study of nonlinear dynamics - the sym-bolic dynamics.Symbolic dynamics originates from the abstract topolog-ical theory of dynamical systems [23], and is gradually de-veloped and perfected in the study of one-dimensional maps[24]. In symbolic dynamics, the state is expressed as an infi-nite sequence of finite set of symbols and the dynamics is rep-resented by symbol shifts, although certain metric informationis lost in this translation [25]. To establish symbolic dynam-ics in a nonlinear system, one of the most critical problemsis how to find a reasonable symbolic partition [26]. In eachpartitioned area, there should be no folding which is typicalin the phase space of a chaotic mapping, where some mono-tonicity or linearity could be envisioned. If the unstable di-rection is one-dimensional, then the partition is realized by aset of "boundary points", determined usually by checking theintersections of stable and unstable manifolds [27–29]. Inter-estingly, we find in this paper that properly selected eigen-functions of the Koopman operator may be used to divide thephase space in a way that is consistent with the symbolic par-tition. This observation is successfully confirmed in several1-d maps with different characteristics as well as in the well-known Hénon map when it is chaotic.The rest of the paper is organized as follows. The con-cept and some properties of the Koopman operator are givenin Sec. II, together with a discussion of phase space partitionbased on its eigenfunctions. Several chaotic maps in one ortwo dimensional phase space are used as examples to demon-strate the implementation of the theory in Sec. III. The cor-respondence between the extrema of eigenfunctions and theboundary points of symbolic partition is checked in great de-tail. In Sec. IV, the whole paper is summarized and possiblefuture directions are pointed out. II. KOOPMAN OPERATOR AND PARTITION OF THEPHASE SPACEA. Koopman Operator and Dynamical Mode
A discrete-time dynamical system, for any point x p ∈ P inthe phase space P , defines x p + = T ( x p ) , where p is the num-ber of iterations, and T gives the iteration law. The Koopmanoperator U acts on a phase space function f ( x ) giving U f ( x ) = f ( T ( x )) (1)which is a linear operator that can be spectrally decomposed.The eigenvalue λ and the eigenfunction φ ( x ) of the Koopmanoperator satisfies U φ ( x ) = φ ( T ( x )) = λφ ( x ) . (2)For example, a one-dimensional discrete dynamical sys-tem x p + = x p gives T ( x ) = x . For an observable function f ( x ) = x , U f ( x ) = f ( T ( x )) = f ( x ) = x = ˜ f ( x ) . That is,the Koopman operator acts on the function and induces thefollowing mappings x U −→ x (3a) f ( x ) U −→ ˜ f ( x ) . (3b)According to Eq. (2), we may obtain two eigenfunctions φ ( x ) = φ = x , and their corresponding eigenvalues are λ = λ = φ ( x ) φ ( x p ) = U φ ( x p − ) = λφ ( x p − )= λ U φ ( x p − ) = λ φ ( x p − ) · · · = λ n − U φ ( x ) = λ n φ ( x ) . (4) We see that the eigenfunctions of the Koopman operatorare related by integer powers of eigenvalues along an orbit.Therefore, these functions are very special in the study ofphase space dynamics of nonlinear systems, which may becalled nonlinear modes in analogy with the eigen-modes ina linear system. If enough modes are identified, the globaldynamics should become clear. However, which are the im-portant modes and how to find them, what kind of roles theyplay are all important questions in the Koopman analysis.In particular, when λ =
1, the eigenfunction φ ( x p ) = φ ( x p − ) = · · · = φ ( x ) ( p ∈ Z is the discrete time index), so theeigenfunction does not change on an orbit, i.e., φ ( x ) = C ( C isa constant). Therefore, different constant eigenfunctions in-dicate different classes of solutions. When | λ | = λ = e i θ ), the eigenfunction | φ ( x p ) | = | φ ( x p − ) | = · · · = | φ ( x ) | ,the modulus of which does not change with time (while thephase changes), and thus corresponds to a constant of mo-tion while the phase increases linearly and leads to oscillationsalong an orbit.For a complex dynamical system, it is often difficult to de-scribe an orbit in the whole phase space. If we divide the phasespace of the dynamical system into distinct regions resultingfrom different dynamical modes, system evolution may be ex-pressed in terms of itinerary of visiting these regions, thatis, it is possible to carry out the dynamical mode decompo-sition (DMD) of the phase space [30]. The eigenvalues andeigenfunctions of the Koopman operator provide us with away for such division. For eigenfunctions with eigenvalues of | λ | =
1, the points with the same modulus belong to an invari-ant set. The work in [31] explores finite-dimensional linearrepresentations of nonlinear dynamical systems by restrictingthe Koopman operator to an invariant set. If the phase spacecan be divided with these invariant sets, we can realize a blockdescription of dynamical modes in the phase space.The idea of describing the phase space with qualitativelydifferent regions reminds us of symbolic dynamics. In a dis-crete dynamical system, we can name different regions of thephase space with different symbols, for example, using thesymbol "0" and "1" to divide the phase space into two differ-ent regions. Each point corresponds to an infinite sequence ofsymbol "0" or "1" according to its visitation of these two re-gions during map iterations. For a permissible finite sequence,a group of points may be identified which have this sequenceas a common visitation itinerary at the initial finite steps. Inthis way, we are able to qualitatively determine how manydifferent orbits are generated in an evolution and estimate theimportance of individual orbits [32]. Even when traditionalanalytical calculations do not apply, we may still try to in-voke a coarse-grained symbolic dynamics to study propertiesof dynamical systems.We call the critical points that divide the phase space as"boundary points". In symbolic dynamics, one of the basic is-sues is how to find "boundary points". As discussed above,specific eigenfunctions of the Koopman operator may helpdistinguish regions with different characteristics, which couldin principle be used to locate these boundary points for a sym-bolic partition.
B. Numerical representation
We obtain the eigenvalues and eigenfunctions of the Koop-man operator by numerical calculation. In order to accu-rately describe an eigenfunction, we can select a group of ba-sis functions g ( x ) , g ( x ) , · · · , g m ( x ) , which make up a finite-dimensional function space, For any function f ( x ) , we mayapproximate it with the expansion f ( x ) = m ∑ i = α i g i ( x ) . (5)Of course, the accuracy of the approximation depends on f ( x ) itself as well as the truncation order m . If we can describethe evolution of all basis functions in this function space, anapproximation of the Koopman operator is obtained.More explicitly, we select the data at the time point p as { x p , x p , · · · , x p n } , and when the time moves to p +
1, thedata evolves to { x p + , x p + , · · · , x p n + } . The basis func-tion g i ( x ) , i = , , · · · , m , as well as the evolved functionsare expressed as column vectors at the known data points { x p , x p , · · · , x p n } , thus forming two data matrices of dimen-sion n × m : K = ( g ( x p ) , g ( x p ) , · · · , g m ( x p ))= g ( x p ) g ( x p ) · · · g m ( x p ) g ( x p ) g ( x p ) · · · g m ( x p ) ... ... . . . ... g ( x p n ) g ( x p n ) · · · g m ( x p n ) , (6a) L = ( g ( x p + ) , g ( x p + ) , · · · , g m ( x p + ))= g ( x p + ) g ( x p + ) · · · g m ( x p + ) g ( x p + ) g ( x p + ) · · · g m ( x p + ) ... ... . . . ... g ( x p n + ) g ( x p n + ) · · · g m ( x p n + ) = ( ˜ g ( x p ) , ˜ g ( x p ) , · · · , ˜ g m ( x p )) . (6b)Each column of K and L is a discrete approximation of a func-tion in the phase space. When the number of points is largeenough, any smooth function with mild oscillation could berepresented with reasonable accuracy. K and L are related bythe Koopman operator K U −→ L . (7)If L is also viewed as a function of x p , the matrix representa-tion of the Koopman operator can be determined by Eq. (7).More explicitly, we choose to write the relation as K ˜ U = L . (8)After obtaining the matrix representation ˜ U of the Koopmanoperator, we can further compute the eigenvalues and eigen-functions.The Koopman operator acts in an infinite dimensional func-tion space. Although our function space dimension is only m,it has been proved [33] that both the spectra measures and pro-jectors of the operators converge to their infinite-dimensionalcounterparts in the limit m → ∞ . After selecting appropriate basis functions and constructingthe data matrices K and L , we obtain m eigenvalues and eigen-functions of the Koopman operator based on Eq. (8). Accord-ing to Eq. (4), we are concerned with the eigenvalue λ =
1, be-cause the corresponding eigenfunction does not change alongan orbit and thus may be used to characterize different ergodiccomponents. However, our calculation is done on a finite pro-jection of the infinite dimensional space, which brings errorsinevitably, so in practice the eigenvalue of λ ≈ C. Function Space: Selection of Basis Functions
According to the discussion in the previous section, to con-struct the data matrix K and L , we have to select a group of ba-sis functions { g i ( x ) } , i = , , · · · , m , which define the functionspace for the approximation and should capture the main dy-namical features of the system under investigation. A properselection will make the computation much easier and allows arobust reconstruction of relevant eigenfunctions, thus consti-tuting an essential part of the whole scheme.For computational efficiency and convenience, usually acomplete set of orthogonal functions are selected. As usual,the orthogonality of the basis functions g i ( x ) is given by h g i , g j i = δ i j , (9)where the inner product h f , h i = R M f ( x ) h ( x ) dx is defined forany two functions f ( x ) and g ( x ) . In a rectangular area, de-pending on the boundary conditions, exponential functions orLegendre polynomials are commonly used orthogonal basisfunctions. However, both of them are non-local functions,which seem hard to capture functions defined on a submani-fold in the phase space. For a dissipative evolution, the asymp-totic attractor often concentrates on a low-dimensional hyper-surface, and a lot of these global functions are needed to depictthis dimension reduction. A more convenient representation isgiven by sets of local functions whose supports are bounded.The following are some examples of basis functions in aninterval on the x − axis g R ( x ) = √ m , ( i − m x < im ) , ( otherwize ) , i = , , · · · , m (10a) g G ( x ) = Cexp − ( x − x i ) d j ! , x i = im − m (10b) g F ( x ) = e ik ( π ) x , k = − m , − ( m − ) , · · · , m − , m (10c) g L ( x ) = r k + P i ( x ) , k = , , · · · , m (10d)where a finite resolution is assumed and signalled by the finiterange of the running indices. Eq. (10a) is a basis of rectangu-lar functions defined in the interval [ , ] , which is orthogonaland local, while the Gaussian basis Eq. (10b) is nearly localand approximately orthogonal, the balance of which is con-trolled by the width d j and the center x i = im − m . For laterreference, the Fourier basis Eq. (10c) in [ , ] and the Legen-dre polynomials Eq. (10d) in [ − , ] are also displayed, whichare orthogonal but non-local.Nevertheless, for high-dimensional systems, all the abovebasis may not be good since the attractor which for finite res-olution may be approximated by an intricately curved hyper-surface needs a lot of these functions to characterize, the num-ber being exponentially increasing with the dimension of thephase space. On the other hand, a heterogeneous distribu-tion of the orbits on the attractor may as well entail a slowconvergence of function expansion and hence requires a largeset of basis functions. One solution to both problems is theutilization of the natural basis constructed from the evolutiondata itself. We may take n points along an orbit to approx-imate a function defined in the relevant region of the phasespace. For example, the n values of the observable x , i.e. , g i ( x ) = ( x i , x i + , · · · , x n + i − ) T defines a discrete representa-tion of a linear function that evolves with time marked withthe index i , according to which Kutz et al [34] proposed analternative view of Koopman (HAVOK) analysis. They putsuch defined functions side by side to form a Hankel matrixas the basis set of functions { g i ( x ) } , i = , · · · , m , whose one-step evolution is obtained by applying the Koopman operatorto g i ( x ) : Ug i ( x ) = U ( x i , x i + , · · · , x n + i − ) T = ( x i + , x i + , · · · , x n + i ) T = g i + ( x ) , (11)which then provides the data matrix K and L in Eq. (6). Asin Eq. (7), the matrix approximation of the Koopman opera-tor is thus obtained. All the basis functions are chosen fromthe given evolution data, which are mainly supported on therelevant part of the phase space and already contains intrinsicproperties of the nonlinear system. Nevertheless, if transientdynamics is also of our concern, the basis functions listed inEq. (10) should be used, which usually involves stable mani-folds in the phase space, as shown in the 2-d example below. III. APPLICATION TO SEVERAL TYPICAL EXAMPLES
The scheme of previous section will be applied to severaltypical chaotic maps including variants of the tent map, thelogistic map in one dimension and the Hénon map in two-dimensional phase space. The eigenvalues and eigenfunctionsof the Koopman operator are calculated using discrete evolu-tion data, which may be checked and used for the partition ofthe phase space.
A. One-dimensional map: the tent and the logistic map
The tent map is a one-dimensional piecewise linear mapdefined on interval [ , ] , named after its functional image re- sembling a tent x n + = f ( x n ) = − (cid:12)(cid:12)(cid:12)(cid:12) x − (cid:12)(cid:12)(cid:12)(cid:12) = x n , x ∈ [ , ) − x n , x ∈ [ , ] (12)which has two fixed points: x ∗ = x ∗ = . The tentmap is a chaotic map with a uniform invariant measure, whichstretches all local segments uniformly and then folds in a sym-metric way. On the other hand, the logistic map x n + = f ( x n ) = γ x n ( − x n ) , x n ∈ [ , ] (13)stretches non-uniformly and the invariant measure is singular.It has two fixed points: x ∗ = x ∗ = − γ as well. Inthe following analysis, we take γ =
1. The Eigenfunctions of Koopman Operator and the Partition ofPhase Space
In the tent and logistic map, we employ the Gaussian ba-sis and take n = m = , , ,
16. Usually we get many eigenvalues and eigen-functions, the number of which depends on m . As we men-tioned in the previous section, we are more concerned withthe eigenvalues close to 1 and the corresponding eigenfunc-tions. We plot the selected eigenfunctions (unless otherwisespecified, when the eigenfunction is complex, we use the realpart) in Fig. 3.In Fig. 3, we choose the eigenvalue of the modulus closestto 1. Each eigenfunction is composed of several wave packets.When m increases, the dimension of our function space alsorises, and the approximation to the Koopman operator gainsaccuracy.Similarly, we can use Eq. (8) to calculate the eigenvaluesand eigenfunctions of the Koopman operator by constructingnatural basis functions. We take n = m = , , , T and its iterations T N : every time the number N increases by one, the function graphs of T N doubles the num-ber of undulations by the process of stretching and folding.The eigenfunctions appear to capture this by wiggling morewhen the number of basis increases.So what do these eigenvalues and eigenfunctions represent?Before discussing this issue, we take the logistic map as anexample to analyze the dynamical process and the "boundarypoints" in symbolic partition.The action of the logistic map can be viewed as a processof stretching and refolding on the interval [0,1], as shown inFig. 1, which consists of two steps: the first is to stretch the AB interval to A ′ B ′ with the points 0,0 . A ′ B ′ interval to A ′′ C ′′ , and maps the points FIG. 1. The dynamical process of the logistic map Eq. (13) is de-composed into two mappings. x n x n + '00' '01' '11' '10' FIG. 2. Symbolic partition of the logistic map Eq. (13). x = . x = x = as a "critical point", and denote the areas as "0" and "1"respectively to the left or the right of the critical point. For anypoint x ∈ [ , ] in the "0" region, its image under the map is apoint that lies either in region "0" or in region "1", according towhich we can further divide the "0" region into "00" and "01" TABLE I. Comparison of extremum points ( m = , , ,
5) andboundary points for different level l as plotted in Fig. 4. m l tent map Eq. (12) logistic map Eq. (13)extremum boundary extremum boundary2 0 0.5000 0.5000 0.5585 0.50003 0 0.5000 0.5000 0.4868 0.50001 0.2502 0.2500 0.1847 0.14640.7500 0.7500 0.8299 0.85364 0 0.5000 0.5000 0.4972 0.50001 0.2503 0.2500 0.1378 0.14640.7500 0.7500 0.8602 0.85362 0.1305 0.1250 0.0518 0.03810.3750 0.3750 0.2938 0.30870.6191 0.6250 0.7107 0.69130.8749 0.8750 0.9485 0.96195 0 0.5000 0.5000 0.4995 0.50001 0.2500 0.2500 0.1472 0.14640.7443 0.7500 0.8522 0.85362 0.1195 0.1250 0.0402 0.03810.3787 0.3750 0.3024 0.30870.6278 0.6250 0.6943 0.69130.8750 0.8750 0.9594 0.96193 0.0597 0.0625 0.0087 0.00960.1894 0.1875 0.0818 0.08420.3033 0.3125 0.2226 0.22220.4268 0.4375 0.4055 0.40240.5625 0.5625 0.5942 0.59750.6861 0.6875 0.7764 0.77780.8172 0.8125 0.9182 0.91570.9560 0.9375 0.9912 0.9904 subregions. Similarly, region "1" can be divided into "11" and"10" subregions. If multiple iterations are considered, we canmark the phase space area with arbitrary precision all the waydown to individual point in the limit of infinite sequence ofsymbols.In the previous discussion, a key problem is how to find"boundary points". In the logistic map, x = is a "boundarypoint". From its preimages, we locate two boundary points forthe next level x = ± √ , which further divides the "0" and"1" region into subregions. From the portraits of the eigen-function in Fig. 3 and Fig. 4, we see that the extremum pointsin the profiles are lying close to the "boundary points" of thesymbolic partition. This is the first evidence that properly se-lected eigenfunctions of the Koopman operator may mirrorthese boundary points.The boundary points of the tent map and logistic map canbe obtained by reverse iteration of the dynamical system, andfor both maps there are two preimages for each point. Wecalculate the boundary points level by level, and list them ac-cording to their levels in Fig. 3. The properly selected eigen-functions of the Koopman operator are plotted with the Gaus-sian basis in Fig. 3 and with the natural basis in Fig. 4. Toguide the vision, we mark the extrema of the eigenfunctionsand the boundary points in these plots, and further list theircoordinates in Table I.For different numbers of Gaussian basis functions as shownin Fig. 3, the local extremum points of eigenfunctions in- x -0.9-0.8-0.7-0.6-0.5-0.4 m=2, x -0.7-0.6-0.5-0.4-0.3 m=4, x m=8, x m=16, (a) tent map Eq. (12) x -0.9-0.8-0.7-0.6-0.5-0.4 m=2, x m=4, x -0.5-0.45-0.4-0.35-0.3-0.25 m=8, x -0.4-0.35-0.3-0.25-0.2-0.15 m=16, (b) logistic map Eq. (13)FIG. 3. Comparison of extremum points (empty squares) and boundary points (colored points, indicating different levels): eigenvalues λ (withmoduli closest to 1, expressed in modulus and angle) and eigenfunctions (black lines) of the Koopman operator ( n = , m = , , , x -1.2-1-0.8-0.6-0.4-0.20 m=2, x m=3, x m=4, x -0.4-0.200.20.40.60.8 m=5, (a) tent map Eq. (12) x -1.2-1-0.8-0.6-0.4-0.20 m=2, x m=3, x m=4, x -0.6-0.4-0.200.20.4 m=5, (b) logistic map Eq. (13)FIG. 4. Comparison of extremum points (empty squares) and boundary points (colored points, indicating different levels): eigenvalues λ (with moduli closest to 1, expressed in modulus and angle) and eigenfunctions (black lines) of the Koopman operator ( n = , m = , , , creases with the number of basis functions. When the num-ber of basis functions is small, the approximation is rough,so the extremum points of eigenfunctions are not in goodagreement with boundary points. However, when we increasethe number of basis functions, we find that the extremumpoints of eigenfunctions become increasingly consistent with the boundary points. This is because the added basis functionsprovide more dynamical undulation details to the eigenfunc-tions, which then more accurately depict the boundary points.For different numbers of natural basis functions as shown inFig. 4, we see that the number of extremum points of eigen-functions is doubled when the number of basis functions in-creases by 1. By comparing the coordinates in Table I, wefind that the newly emerged extremum points correspond tothe boundary points at the next level, and with the increase ofthe basis, the discrepancy between extremum and boundarypoints also decreases.We have thus preliminarily verified the consistency be-tween the extremum points of eigenfunctions and the bound-ary points of the symbolic partition. In the previous discus-sion, we learned that these boundary points are actually thecritical points of symbolic dynamics. We believe that this ruleapplies to other cases as long as the chaotic dynamics is origi-nated from local stretching and global folding. As we increasethe number of basis functions, more partition points at differ-ent levels may be identified with greater accuracy, thus en-abling us to have a better understanding of the characteristicsof the dynamics.
2. Robustness of the partition in the presence of noise
In the real world, the evolution of a nonlinear system isoften contaminated by the omnipresent noise. Heninger pro-posed the finest state-space resolution that can be achieved ina physical dynamical system is limited by noise [35], whichin fact revives a form of perturbation theory in [36]. The fi-nite approximation of the Koopman operator presented aboveof course brings error in the computation but at the same timegrants robustness even in the presence of noise. We will inves-tigate what kind of role noise plays and how much the spectralproperties of the Koopman operator are altered. For simplic-ity, the noise is assumed to be Gaussian white. The param-eters in the following computation are: the number of points n = m = , , , σ = . L from the same K is used for the constructionof Eq. (8). Similar procedure is used for the computation asbefore, and we take the eigenvalues closest to 1 to plot thecorresponding eigenfunctions in Fig. 5.Comparing Fig. 3 with Fig. 5, we find that the set of bound-ary points by Koopman operator is almost the same, whichshows the robustness of our scheme.
3. Consistency in different truncations
With the increase of the number of basis functions, theeigenfunctions describe the boundary points more finely: themore the basis functions, the higher the level and finenessof the boundary points prescribed by the eigenfunctions. Inorder to verify our conclusion, we make the following com-parison: when the number of basis functions doubles, we ex-plore the relationship between the extremum points of per-tinent eigenfunctions. In order to compare the similarity indynamical characteristics, we define a correlation coefficientof the eigenfunctions ρ ( φ , φ ) = ∑ i ( φ i − φ )( φ i − φ ) q ( ∑ i ( φ i − φ ) )( ∑ i ( φ i − φ ) ) , (14) where φ , φ are the corresponding eigenfunctions, and φ , φ are their mean values.To ensure the consistency of corresponding eigenfunctionsfor different truncations, the following scheme is designed tosearch φ m when φ m is knownarg min φ m || U φ m − λφ m || + µ | ρ ( φ m , φ m ) | , (15)where U is the Koopman operator and λ is the eigenvalueof φ m . The first term in Eq. (15) is the definition of eigen-function and the second term is the correlation coefficient de-fined in Eq. (14). µ is a relative weight. Generally speaking,we should ensure first the correctness of the eigenfunction, sousually µ ≪ m = m =
16 and inFig. 6 compared the results (blue curves for m = m =
4. One-dimensional maps with multiple peaks
We have successfully partitioned the phase space of 1-dunimodal maps with eigenfunctions of the Koopman operator,but its applicability to other systems needs further investiga-tion. Here we discuss certain multimodal maps, as given inEq. (16a) and Eq. (16b). f ( x ) = x , ( x . ) − x , ( . x . ) x − , ( . x . ) − x , ( . x ) (16a) f ( x ) = x , ( x . ) − x , ( . x . ) x − , ( . x . ) . − x , ( . x ) (16b)Fig. 7 displays the relevant eigenfunctions of the Koopmanoperator for f and f . In both maps, there are three criticalpoints: 0 .
25, 0 . .
75 at the first level. The positions ofthe boundary points seem to coincide quite well with those ofthe extremum points of the corresponding eigenfunction.However, there are some differences in the two maps. Forexample, in map f , the function values at x = .
25 and x = .
75 are twice as different, but the eigenfunctions of thetwo maps are surprisingly similar, which had made us doubtthe accuracy of our computation. By checking more eigen-function plots, we find that this is true only for the principle x -0.9-0.8-0.7-0.6-0.5-0.4 m=2, x m=4, x -0.5-0.45-0.4-0.35-0.3-0.25 m=8, x m=16, (a) tent map Eq. (12) x -0.9-0.8-0.7-0.6-0.5-0.4 m=2, x m=4, x -0.5-0.45-0.4-0.35-0.3-0.25 m=8, x m=16, (b) logistic map Eq. (13)FIG. 5. Same as in Fig. 3, but subject to Gaussian white noise with variance σ = . x m=8m=16 (a) tent map Eq. (12) x -0.038-0.036-0.034-0.032-0.03-0.028-0.026-0.024-0.022-0.02-0.018 m=8m=16 (b) logistic map Eq. (13)FIG. 6. Correspondence between eigenfunctions with different m ( n = m =
16 (red) starting from itscounterpart with m = µ = . eigenfunction (with eigenvalue closest to 1). If we plot thesecond most relevant function, as in Fig. 8, the two eigenfunc-tion profiles are very different. This shows that other eigen-functions can describe dynamical difference between the two,and the eigenfunctions for disparate eigenvalues contain dif-ferent information of the dynamics. B. Two-dimensional map: Hénon map
Hénon map is a discrete mapping system generating chaoswith appropriate parameters [37], which reads x n + = y n + − ax n y n + = bx n (17)and has two fixed points: x ∗ = b − ± p ( b − ) + a a , y ∗ = bx ∗ . For the parameter value a = . b = .
3, Hénon map x -0.9-0.8-0.7-0.6-0.5-0.4 m=2, x -0.7-0.6-0.5-0.4-0.3 m=4, x m=8, x m=16, (a) f map Eq. (16a) x m=2, x -0.7-0.6-0.5-0.4-0.3 m=4, x -0.5-0.45-0.4-0.35-0.3-0.25 m=8, x -0.35-0.3-0.25-0.2-0.15 m=16, (b) f map Eq. (16b)FIG. 7. Comparison of extremum points (empty squares) and boundary points (colored points, indicating different levels): eigenvalues λ (for the first eigenvalue, expressed in modulus and angle) and eigenfunctions (black lines) of the Koopman operator ( n = , m = , , , f map Eq. (16a) and (b) f map Eq. (16b). x -0.6-0.4-0.200.20.4 m=2, x -0.6-0.4-0.200.20.40.6 m=4, x -0.4-0.200.20.40.6 m=8, x -0.3-0.2-0.100.10.20.3 m=16, (a) f map Eq. (16a) x -0.500.51 m=2, x -0.6-0.4-0.200.20.4 m=4, x -0.4-0.200.20.40.60.8 m=8, x -0.4-0.200.20.4 m=16, (b) f map Eq. (16b)FIG. 8. Comparison of extremum points (empty squares) and boundary points (colored points, indicating different levels): eigenvalues λ (forthe second eigenvalue, expressed in modulus and angle) and eigenfunctions (black lines) of the Koopman operator ( n = , m = , , , f map Eq. (16a) and (b) f map Eq. (16b). produces chaos and a strange attractor appears in the phasespace.In a two-dimensional map, the calculation of the Koop-man operator eigenfunction is slightly different from the one-dimensional counterpart. For example, we may use two- dimensional Gaussian functions as basis g ( x , y ) = Cexp − ( x − x j ) − ( y − y j ) d j ! , (18)where C is a normalization constant, ( x j , y j ) is the center ofthe Gaussian wave packet, and d j is the width of the wave0packet. In the rectangular area [ − . , . ] × [ − . , . ] whichcontains the attractor, we deploy the Gaussian basis functionswith centers uniformly distributed on a 50 ×
50 lattice and awave packet width d j = /
45, and with the number of points n = ×
100 calculate the eigenvalues and eigenfunctions.We take the first four eigenfunctions whose eigenvalues areclosest to 1. When the eigenfunction value is complex, thereal part is used for the plot. Fig. 9 displays the color imagesof the eigenfunctions. Note that the profile of the eigenfunc-tion is very similar to that of the stable manifold of Hénonmap [29]. The boundary points of the Hénon map could bechosen from the homoclinic tangencies of the stable and un-stable manifold [28, 29]. Here, the asymptotic state of a typi-cal point is on the closure of the unstable manifold while theKoopman eigenfunctions are capable of revealing the stablemanifold, which partly explains why the partition boundariescould be located in this new approach.Alternatively the Koopman operator may be approximatedwith the natural function basis. However, it is slightly differ-ent from the one-dimensional case and involves a minor com-plication. In the two-dimensional map, there are two state ob-servables x and y , either of which can be used for the construc-tion. For finite resolution, different observables may empha-size different aspects of the dynamics and thus have differentprojections on the eigenfunctions. In the discussion of Hénonmap, we take the x − direction as an example to emphasize thestretching and folding dynamics. At the same time, for thisset of parameter values, the Hénon map is dissipative and thestrange attractor has a fractal dimension slightly greater than1. The natural basis is just defined on the attractor and weonly plot the eigenfunctions supported on the attractor.Fig. 11 shows the eigenfunction plots on the attractor with n = m = , , ,
4. Since the attractor is thin and dif-ficult to represent, we use a three-dimensional plot, portrayingthe eigenfunction of the Koopman operator with height andcolor at the same time, and also project it to the x − y plane.As can be seen from Fig. 11, with the increase of the num-ber of basis functions, the structure of eigenfunction imagesbecomes increasingly complicated, and the number of ex-tremum points increases exponentially, which is consistentwith our conclusion in one-dimensional mapping. For sim-plicity in a coarse-grained picture, we regard the attractor ofHénon map as consisting of four main segments over whicheigenfunctions are plotted. When m =
1, the eigenfunctionidentifies only the principle group of critical points, one oneach main segment and dividing the attractor of the Hénonmap into two parts. When m =
2, there is an extra set of ex-tremum points, which divide each subsegment from the previ-ous step into two parts, which is very similar to the 1-d case.The numerical values of boundary points are given in [28],the principle four of which is listed in Table II. They are alllocated in the main bending region of the attractor, and closeto the extremum points identified in the eigenfunction. Inorder to find the boundary points at other levels, we carryout backward evolution of the four boundary points to obtaintheir preimages. As depicted in Fig. 10, we label the preim-ages with positive numbers, and 0 represents the original fourboundary points.
TABLE II. Comparison of extremum points ( m = , , ,
4) andboundary points (given in [28]) at different level l of the Hénon map( a = . , b = .
3) as plotted in Fig. 11. m l extremum boundary1 0 ( 0.7079, 0.0120) ( 0.7021,-0.0044)( 0.8019, 0.0191) ( 0.7986, 0.0019)( 1.2326,-0.0153) ( 1.2307,-0.0249)( 1.2729,-0.0091) ( 1.2717,-0.0205)2 0 ( 0.7052,-0.0009) ( 0.7021,-0.0044)( 0.7995, 0.0041) ( 0.7986, 0.0019)( 1.2720,-0.0192) ( 1.2307,-0.0249)( 1.2720,-0.0192) ( 1.2717,-0.0205)1 (-0.1418,-0.3100) (-0.0148,-0.2976)(-0.1258,-0.2206) ( 0.0065,-0.2013)(-0.2164, 0.2936) (-0.0832, 0.2403)(-0.2164, 0.2936) (-0.0682, 0.2782)3 0 ( 0.7045,-0.0026) ( 0.7021,-0.0044)( 0.7995, 0.0041) ( 0.7986, 0.0019)( 1.2713,-0.0217) ( 1.2307,-0.0249)( 1.2713,-0.0217) ( 1.2717,-0.0205)1 (-0.0330,-0.2981) (-0.0148,-0.2976)(-0.0121,-0.2041) ( 0.0065,-0.2013)(-0.0982, 0.2431) (-0.0832, 0.2403)(-0.0862, 0.2802) (-0.0682, 0.2782)2 (-1.0459, 0.3549) (-0.9918, 0.3625)(-0.7503,-0.3700) (-0.6711,-0.3630)( 0.9130,-0.1574) ( 0.8011,-0.1846)( 1.0081, 0.1176) ( 0.9275, 0.1361)4 0 ( 0.7045,-0.0026) ( 0.7021,-0.0044)( 0.8002, 0.0059) ( 0.7986, 0.0019)( 1.2713,-0.0217) ( 1.2307,-0.0249)( 1.2713,-0.0217) ( 1.2717,-0.0205)1 (-0.0086,-0.2954) (-0.0148,-0.2976)( 0.0138,-0.2002) ( 0.0065,-0.2013)(-0.0856, 0.2413) (-0.0832, 0.2403)(-0.0724, 0.2787) (-0.0682, 0.2782)2 (-0.9937, 0.3495) (-0.9918, 0.3625)(-0.6805,-0.3638) (-0.6711,-0.3630)( 0.8132,-0.1780) ( 0.8011,-0.1846)( 0.9338, 0.1347) ( 0.9275, 0.1361)3 ( 1.2205, 0.0336) ( 1.2084, 0.0523)(-1.2334, 0.3797) (-1.2099, 0.3783)(-0.5413, 0.3240) (-0.6154, 0.3313)( 0.3897, 0.2237) ( 0.4535, 0.2154)
The eigenfunctions together with the boundary points areplotted in Fig. 11. The coordinates of extremum points ofeigenfunctions and boundary points are listed in Table II. Bychecking their positions in Fig. 11 and comparing their coor-dinates in Table II, we find that the extremum points of eigen-functions always correspond to the boundary points of thesymbolic partition of the Hénon map, and higher-level bound-ary points emerge with the increase of the number of basisfunctions. As such, the extremum points of the eigenfunctionmatch better with the corresponding boundary points, whichis consistent with the one-dimensional case.1
FIG. 9. Eigenvalues λ (the first four eigenvalues with moduli closest to 1, expressed in modulus and angle) and eigenfunctions (color plots) ofHénon map Eq. (17) ( a = . , b = .
3) in the rectangular area [ . × . ] with the Gaussian basis ( n = , m = , d j = / IV. CONCLUSION
Koopman operator describes the evolution of observablefunctions in phase space, which could be used to excavateimportant patterns of dynamics. In a spectral decompositionof the Koopman operator, numerous eigenvalues and eigen-functions are identified but it remains a big challenge to iden-tify the most relevant ones and reveal their dynamical signif-icance. In the current paper, we find that extrema of properlyselected eigenfunctions overlap with the boundary points ofsymbolic partitions in the 1-d or 2-d maps with chaotic dy-namics. With well-chosen basis functions, the spectral prop-erties of the Koopman operator could easily be computed andthe accuracy increases with the truncation order. Comparedto other basis, the natural one directly drawn from the evolu-tion data is the most efficient, which should be applicable to a wide-variety of other non-linear systems.Symbolic partition is always a mind-boggling problem inthe study of nonlinear dynamics, especially when the phasespace is high dimensional [38–40]. Most previous approachesare concentrated in the construction of stable and unstablemanifolds or identification of some topological index, whichare hard to invoke when orbit structures are complicated, be-ing usually the case for coupled nonlinear systems. In the cur-rent paper, a totally different route is taken by considering thespectral properties of the Koopman operator, thus avoiding thenearly impossible description of possibly intricate geometricfeatures in the phase space. Therefore, it is reasonable and de-sirable to extend our new scheme to the treatment of complexnonlinear systems and achieve proper symbolic description.Although the analysis is demonstrated in the well-know 1-d or 2-d maps, further investigation is still needed for its ex-2 -1.5 -1 -0.5 0 0.5 1 1.5 x -0.4-0.3-0.2-0.100.10.20.30.4 y FIG. 10. Boundary points of Hénon map Eq. (17): the principle ones(marked by 0) and their preimages (marked by positive integers fordifferent levels). tension, especially in the selection of proper eigenfunctions.In this paper, eigenfunctions with eigenvalues close to 1 arechosen for the symbolic partition. It would be interesting tocheck the roles of other eigenfunctions and propose a sys-tematic scheme for the selection of relevant eigenfunctions.We also checked and compared eigenfunctions for differenttruncations in the computation and found that the extremum points are consistent in these eigenfunctions but new onesemerge with an increasing number of basis functions, unfold-ing higher levels of partition boundaries. Therefore, it seemspossible to utilize the truncation order to coarsen or refine thedescription when necessary, which should be explored in moredetail in future.In a 1-d example in the above discussion, small Gaussianwhite noise is added to the evolution but the partition based onthe Koopman operator seems intact, which shows the robust-ness of the current scheme. It would be interesting to checkif and how the partition changes with the noise intensity sincethere should exist an optimal partition for noise with a finitesize [35, 36]. In the above 2-d example, the x − coordinate isused for the Koopman analysis since the stretching and fold-ing mechanism is most clear in this direction. In a generalsituation, we have no idea whatsoever which are importantdirections. It would be good if a systematic procedure is de-veloped for this purpose. In practice, very often, only partialdata is collected for a system and the desired one may not evenbe included. The problem of how much and what kind of in-formation could be extracted from such data seems to be ofpractical importance. ACKNOWLEDGMENTS
This work was supported by the National Natural ScienceFoundation of China under Grants No. 11775035 and No.11375093, and also by the Fundamental Research Funds forthe Central Universities with Contract No. 2019XD-A10. [1] M. H. Hassoun et al. , Fundamentals of artificial neural net-works (MIT press, 1995).[2] S. Strogatz,
Nonlinear dynamics and chaos: with applicationsto physics, biology, chemistry, and engineering (studies in non-linearity) (CRC Press, 2001).[3] P. Holmes, J. L. Lumley, and G. Berkooz,
Turbulence, Coher-ent Structures, Dynamical Systems and Symmetry. (CambridgeUniversity Press, 1996).[4] L. Landau,
Statistical Mechanics (Butterworth-Heinemann,1973).[5] I. Prigogine,
Non-equilibrium statistical mechanics (CourierDover Publications, 2017).[6] B. O. Koopman, Proc. Natl. Acad. Sci. U.S.A. , 315 (1931).[7] P. Gaspard, G. Nicolis, A. Provata, and S. Tasaki, Phys. Rev. E , 74 (1995).[8] P. Gaspard and S. Tasaki, Phys. Rev. E , 056232 (2001).[9] M. Budiši´c, R. Mohr, and I. Mezi´c, Chaos , 047510 (2012).[10] I. Mezi´c and A. Banaszuk, Physica D , 101 (2004).[11] I. Mezi´c, Nonl. Dyn. , 309 (2005).[12] Y. Susuki and I. Mezi´c, IEEE Trans. Power Syst. , 1894(2011).[13] Y. Susuki and I. Mezi´c, IEEE Trans. Power Syst. , 1182(2012).[14] B. Eisenhower, T. Maile, M. Fischer, and I. Mezi´c, Proc. Sim-Build , 434 (2010).[15] M. Georgescu, B. Eisenhower, and I. Mezi´c, Proc. SimBuild , 40 (2012).[16] P. J. Schmid, L. Li, M. P. Juniper, and O. Pust, Theor. Comput.Fluid Dyn. , 249 (2011).[17] S. Bagheri, J. Fluid Mech. , 596 (2013).[18] I. Mezi´c, Annu. Rev. Fluid Mech. , 357 (2013).[19] P. J. Schmid, J. Fluid Mech. , 5 (2010).[20] M. Korda, M. Putinar, and I. Mezi´c, Appl. Comput. Harmon.Anal. , 599 (2020).[21] H. Arbabi and I. Mezic, SIAM J. Appl. Dyn. Syst. , 2096(2017).[22] Y. Lan and I. Mezi´c, Physica D , 42 (2013).[23] M. Morse and G. A. Hedlund, Am. J. Math. , 815 (1938).[24] J. Crutchfield and N. Packard, Int. J. Theor. Phys. , 433(1982).[25] C. Robinson, Dynamical systems: stability, symbolic dynamics,and chaos (CRC press, 1998).[26] B.-l. Hao, Physica D , 161 (1991).[27] O. Biham and W. Wenzel, Phys. Rev. Lett. , 819 (1989).[28] L. Jaeger and H. Kantz, J. Phys. A: Math. Gen. , L567 (1997).[29] P. Grassberger and H. Kantz, Phys. Lett. A , 235 (1985).[30] A. Alla and J. N. Kutz, SIAM J. Sci. Comput. , B778 (2017).[31] S. L. Brunton, B. W. Brunton, J. L. Proctor, and J. N. Kutz,PLoS One (2016).[32] P. Cvitanovi´c, R. Artuso, R. Mainieri, G. Tanner, and G. Vattay, Chaos: Classical and Quantum (ChaosBook.org, 2020).[33] N. Govindarajan, R. Mohr, S. Chandrasekaran, and I. Mezic, FIG. 11. Comparison of extremum points (black dots) and boundary points (colored dots, indicating different levels): eigenvalues λ (withmoduli closest to 1, expressed in modulus and angle) and eigenfunctions (colored lines) of the Koopman operator ( n = , m = , , , a = . , b = . , 1454 (2019).[34] S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, and J. N.Kutz, Nat. Commun. , 1 (2017).[35] J. M. Heninger, D. Lippolis, and P. Cvitanovi´c, Phys. Rev. E ,062922 (2015).[36] J. M. Heninger, D. Lippolis, and P. Cvitanovi´c, Commun. Nonl. Sci. Numer. Simul. , 16 (2018).[37] C. Simó, J. Stat. Phys. , 465 (1979).[38] A. Gonchenko and S. Gonchenko, Physica D , 43 (2016).[39] A. Gonchenko, S. Gonchenko, A. Kazakov, and D. Turaev, Int.J. Bifurc. Chaos , 1440005 (2014).[40] S. V. Gonchenko, I. Ovsyannikov, C. Simó, and D. Turaev, Int.J. Bifurc. Chaos15