Ordinal Poincaré sections: Reconstructing the first return map from an ordinal segmentation of time series

The Poincaré map (also known as the first return map) is a useful tool for studying dynamical systems—projecting a continuous dynamical system to a lower dimensional map. A Poincaré section is the intersection of an orbit in the state space of a continuous dynamical system with a particular lower-dimensional surface transverse to the system’s flow.¹ More specifically, consider an orbit with initial conditions within a section of space, which then leaves that section, and observe the point at which this orbit first returns to the section. The first return map (FRM) is then created by making a map that connects the first point in the section to the second, the second point to the third, and so on. The transversality of the Poincaré section means that the orbits begin in the subspace and flow through it rather than parallel to it. Hereafter, we will refer to the $n−1$ dimensional surface transverse to the flow as the Poincaré section, and the map we construct from successful intersections with that section as the First Return Map.

Following the idea of the Poincaré recurrence theorem, Kakutani introduced the FRM in 1943 with the goal of investigating the properties of induced measure-preserving transformations.² A Poincaré map is a discrete dynamical system with dimension $n−1$ derived from the original continuous dynamical system with dimension $n$⁠. It is frequently used to simplify the analysis of the dynamical system because it preserves many of the original system’s properties by transforming periodic and quasiperiodic orbits into periodic and quasiperiodic points.

These FRMs are useful for studying oscillatory dynamics and have been used to explain the origins of chaotic dynamics in specific models.³  Figure 1(a) depicts the surfaces of maxima of $x$ with $x ˙=0$ in gray as a Poincaré section on the Lorenz attractor. The attractor’s entry points to this section are marked with dots. One of the $x ˙=0$ surfaces crosses the attractor in a negative $x$ region (red dots), while the other passes through a positive $x$ region (black dots). Figure 1(b) shows a portion of the Lorenz time series after passing through the transient. Red and black dots highlight the local maxima entering the $x ˙=0$ Poincaré section. Figure 1(c) depicts the FRM generated by red and black dots. The concept of constructing the FRM from the points entering a Poincaré section was first introduced in a similar figure in 1963 for the Lorenz attractor’s $z$ time series in Lorenz’s famous paper.⁴ 

Despite all of the desirable features of the first return map, it is not always easy to find an ideal Poincaré section that can accurately capture the system’s characteristics.⁵ Also, suppose the system’s equations are not fully available, and only its output time series is known. In that case, in order to use the Poincaré section, the time series must be transformed into a suitable higher-dimensional structure that allows the application of standard high-dimensional analysis methods.⁶ One of the most popular methods to do so is the reconstruction of phase space from time series in a way that is equivalent to a system’s original phase space.⁷ Takens’ delay embedding theorem for an infinite, noise-free time series to unfold the attractor of its generator system is a paradigmatic example of embedding techniques.⁸ This method is widely applied in practice; however, it is not always easy to find the appropriate dimension and lag in the embedding solution that produces an attractor topologically similar to the original system and is computationally efficient.⁹ 

In addition to embedding, recent approaches to convert time series data to a complex network form have gained popularity and allow one to examine time series in higher dimensions and use complex network analysis tools to extract time series features.^10–12 Several algorithms can generate a complex network representation from a time series.^12–15 The one that we are going to use in this study is the ordinal partition network approach.¹⁶ Converting a continuous time series (albeit sampled at discrete times) to a complex network involves two important approximations. Nodes of the networks are associated with regions of state space, and so there is coarse-grained quantization of original observations. Furthermore, the network representation encodes probability-based transitions between nodes. It allows one to compute properties related to the underlying invariant measure but does not preserve the deterministic dynamics.¹⁷ 

The ordinal partition network (OPN) was originally introduced in 2013.¹⁶ This dynamical symbolic mapping¹⁷ enables the approximation of deterministic dynamics by examining the statistical properties of a stochastic model. In simpler terms, this partitioning defines a map between the time delay embedded of the observed data and the symbolic dynamics of interest. This method is computationally efficient since for a given choice of parameters, only the weights of the links increase and no additional symbols are added for long time series.¹⁸ 

Unlike OPN approaches, which divide state space into distinct regions and provide a stochastic approximation to the dynamics, the ordinal-based FRM approach that we propose here preserves the original fidelity while constructing the deterministic discrete map. We construct it by assigning an ordinal partition to each point in the time series and taking the change in the ordinal partition between two nearby points as the entrance to that specific ordinal partition section.

In our paper, we present a generic method for detecting dynamical behavior of a time series by treating ordinal partitions as Poincaré sections. Although these sections differ slightly from the primary definition of Poincaré sections, we show that they can behave similarly and may be used to extract the main dynamical features of the system. This method extracts features from a single time series without requiring access to the entire attractor or performing full embedding. Furthermore, by combining the FRMs from different ordinal partitions, additional information regarding the dynamical features can be elucidated. It is also noise resistant as it is built on ordinal partitions, which allows it to be used on real-world data with significant noise contamination.¹² We used a numerically generated time series from the Lorenz system to validate our theoretical discussion. We assign a symbol to each point of the time series by ordinal partitioning the amplitude of the time series points in each window. For the next step, for each set of ordinal partitions, we create an OPN from all the points in it to calculate the entropy for each set. We can determine which ordinal partitions are good Poincaré sections by ranking them based on their OPN’s entropies and comparing their FRMs to those obtained by the original forms of Poincaré sections.

The rest of this paper is organized as follows. In Sec. II, we first introduce the methods and measures used to encode the time series in a symbolic sequence, which allows one to estimate the dynamical features of the time series. Then we introduce various types of permutation entropies for ranking ordinal partitions. The results of applying this method to the Lorenz time series are presented in Sec. III, and the conclusion is presented in Sec. IV. The results of this method’s investigation on Lorenz time series with other ordinal parameters as well as time series from some other dynamical systems are given in the  Appendix.

Let $x t$ represent a scalar time series measured from a dynamical system with known or unknown equations. We generate an OPN from this time series in the following manner. We consider an ordinal window of length $ℓ$ and slide it over the time series. The ordinal window length $ℓ$ determines the range of time series points to be analyzed. By dividing this window length into $m−1$ equal parts with length $τ$⁠, we have $m$ equally spaced points to compare their amplitudes and assign an ordinal partition to them. In other words, the ordinal window is made up of $m$ points, with $τ$-points gap between each pair, and the window length is $ℓ=(m−1)τ$⁠. We note that all of the parameters $ℓ$⁠, $m$⁠, and $τ$ are integers, and the values of $ℓ$ and $m$ should be chosen in such a way that $τ= ℓ m − 1$ becomes an integer as well. Figure 2 illustrates a simple example for $ℓ=4$⁠, $m=3$⁠, and $τ=2$⁠. A fourth parameter is $w$⁠, which represents the number of non-overlapping points of neighbor windows. In this study, we chose $w=1$⁠, which means that the windows have the most overlap, and we need $N−(m−1)τ$ ordinal windows to analyze the entire $N$-point measurement time series.

In each ordinal window, the corresponding ordinal partition $o ( i )=( π 1, π 2,…, π m)$ where $π j∈{1,2,…,m}$ (⁠ $π j≠ π k if j≠k$⁠) could be defined for each vector $z ( i )= ( x i , x i + τ , … , x i + ( m − 1 ) τ )$⁠. In this way, $o ( i )$ captures the relative order of $z ( i )$ vector elements. This symbolic ordering can be done in two different ways: amplitude ranking or chronological index ranking.¹⁷ Amplitude ranking ranks the points based on the relative amplitude, so the point with the lowest amplitude has rank $m$ and the point with the highest amplitude has rank $1$⁠. On the other hand, chronological index ranking uses the time index to sort the points such that it is the index of the lowest value that is placed first in $o$⁠. Table I shows these two ways of ranking for $m=3$⁠. In this work, we used chronological index ranking for creating the OPN.

Using chronological index ranking in each window, we can sort each window’s selected $m$ points into at most $m!$ different orders as described above. Since some of the possible ordinal rankings may never occur within the time series, the network has at most $m!$ nodes. We can build a network for each time series in the following way. Associating each order of $z ( i )$ to a network node allows for considering a weighted directed graph. The nodes of this network are different occurrences of ordinal partitions of $m$ selected points in each window. So, considering larger $m$ in each window results in more nodes in the final network. The links between nodes illustrate the changes from each state to another. We also allow for self-loops when the state returns to itself. This network demonstrates how ordered patterns in the sampled time series have changed over time. Figure 2 depicts this structure on a sampled time series of a map. The relative magnitude of $m$ points and their indices are taken into account in each window, and the adjacency matrix of the dynamical network is then created by arranging sequential ordinal symbols.

A reasonable prescription for selecting the most appropriate ordinal window length (⁠ $ℓ$⁠) for extracting dynamical features of the time series is to use traditional embedding methods. This is because both tasks (reconstruction of phase space from the output time series, which is equivalent to a system’s original phase space, and creating an OPN from a time series) practically revolve around the same central issue of determining how long the time series in a window can preserve the topological properties of the signal. If $ℓ$ is too small, it will incorrectly observe very small noise fluctuations as system dynamics. If $ℓ$ is too large, it will fail to recognize the fine scale structure of the time series. Any of the standard embedding methods for selecting embedding lags can be used to find a suitable $ℓ$⁠.^19–22 In this work, we use the time delay method provided by Takens⁸ to select the appropriate ordinal window length (⁠ $ℓ$⁠). More details are given in Sec. III.

Takens’ embedding theorem is not the only possible solution, and in some cases, it does not work at all. For example, this method cannot be applied to delay systems with infinite dimensions. Also, in other non-delayed systems, which is challenging to find the appropriate dimension and delay for embedding the system in higher dimensions due to the complexity of their dynamics, any alternative method can be used to find a suitable window length for sequential partitions. Therefore, despite the fact that finding a suitable ordinal window is related to the length of the embedding window, ordinal partitioning has fewer limitations than embedding methods.

Because each set of points in ordinal partitioning is assigned to one and only one symbol, the ordinal partitions divide the system’s attractor into several regions such that the union of all these regions covers the entire attractor. The border of each of these regions is defined by a set of points that separate these areas. We assert that each of these boundaries can be considered a Poincaré section on the attractor, and the FRM can be studied on it. Depending on which part of the attractor this ordinal partition cut is on, it can be considered as an appropriate or inappropriate Poincaré section on the attractor.

To calculate the entropy for each set of ordinal partitions’ points, we select all the points in the initial time series that have the same ordinal partition and combine them to form a new time series. The entropy for each ordinal partition is then calculated by comparing the amplitudes of the points in this new time series with each other and generating an OPN from it. We consider window options $m ′=3$ and $τ ′=1$ for creating the OPN from each ordinal partition’s time series. In this regard, for the first time, we proposed two new types of entropy that connect the probability of occurrence of each set of ordinal partition’s points in the first time series and the OPN created from each of them.

To divide the time series into different ordinal partitions, we consider an ordinal window as shown in Fig. 3 and move it through the entire time series. This allows us to assign a symbol to each time series point. We then generate an OPN from each of these symbols to determine $h w$ and $h w t$ of each collection. Figure 4 depicts the time series of the descending node (ordinal partition [4,3,2,1]) as an example of a new independently constructed time series from a specific ordinal sequence that has been extracted from Lorenz’s time series. To compute $h w$ and $h w t$ for this new time series, we must first create an OPN out of it. The newly constructed time series can be considered as a mapping that samples the original time series and represents the information related to the ordinal sequence in a lower dimension. As a result, the shortest possible window length with no gaps between its points is assumed to build the OPN for calculating the permutation entropy. The ordinal window parameters for constructing the OPN from all sets of ordinal points in the main time series will be $m ′=3$⁠, $τ ′=1$⁠, and $ℓ ′=( m ′−1) τ ′=2$⁠, where we use primes to denote the time series of the ordinal partition of interest. Then, to calculate the scaling factors for each specific ordinal partition, we set $O$ equal to the length of the time series segment and set $O ^$ equal to the number of large points in it that specify the first entries to this section. As previously stated, $T$ and $T ^$ are calculated from the entire time series of the system and are the same for all ordinal partitions, not just the partition under study.

Figure 5(a) displays the $h w$ for each ordinal partition, which is sorted by highest to lowest entropy in the first column. The $h w t$ is displayed similarly in the second column. A color has been assigned to each ordinal partition, which is more clearly shown in the third column. This allows us to assign a color to each time series point, as seen in Fig. 5(b). With this colored time series, we may recreate the attractor in three dimensions according to Takens’ embedding theorem, yielding Fig. 5(c). By comparing the topological properties of this attractor to those of the Lorenz original attractor in Fig. 1(a), it can be seen that the embedded attractor mimics the primary features of Lorenz, such as its two wings and two holes. As can be seen, each ordinal partition has a distinct location on the attractor; thus, the entry points of ordinal partitions with high entropy values can be thought of as Poincaré sections on the attractor. On the other hand, the ordinal partitions with low entropy values are just a few single dots on the time series, randomly arranged on the embedded attractor, and cannot be considered as a good section that captures the dynamics of the system.

Attractor regioning by ordinal partitions occurs in all dynamical systems, regardless of their dimension or complexity. The  Appendix shows similar results for Rössler and Mackey–Glass systems as examples. It is obvious that because Mackey–Glass has a more complicated dynamical behavior, it experiences a richer collection of ordinal partitions.

Figure 6(a) depicts the value of $h w t$ (black line) and $h w$ (green line) of each ordinal partition, which are sorted in the descending order based on $h w t$⁠. The range of $h w t$ is shown on the left axis, and the range of $h w$ is shown on the right axis. Although the ranges of these two entropies differ, both show three levels of entropy. The levels of $h w t$ are shown in solid boxes, and the levels of $h w$ are shown in dashed boxes. In both, blue boxes represent the first level of entropy, yellow boxes represent the second level, and red boxes represent the third level. The graph illustrates that the first and second entropy levels are switched in different entropy methods, but the third level has the same group of ordinal partitions. At the first level, $h w$ has a higher value than $h w t$⁠, but $h w t$ has a higher value at second and third levels. Generally, the variance of $h w$ is greater than the variance of $h w t$⁠. This is because the weight of the ascending and descending ordinal partitions [jade green and light green parts in Fig. 5(c)] in $h w$⁠, which is calculated by the number of all points that have that ordinal pattern, is much higher than the weight of all other ordinal partitions. However, the $h w t$ calculates the weights only based on the entry points to each section, putting them in the same range.

As shown in the  Appendix, the number of levels and their entropy amplitude in both $h w$ and $h w t$ can be changed due to the system’s dynamics and the number of points considered in each window of the ordinal partition. The changes in $h w t$ can be smoother as more points are considered in the ordinal window, allowing the first and second levels of entropy to be merged with each other. Furthermore, depending on the system’s dynamics, the number of levels in $h w$ and $h w t$ may differ, or there may be no obvious levels at all.

Figure 6(b) depicts a similar procedure when noise is incorporated. Random Gaussian noise with varying signal-to-noise ratios (SNR) has been added to the signal. The following SNRs have been considered: 100 (which is nearly equal to no noise), 70, 60, 50, 40, 30, 20, and 10 dB. At each level of noise, corresponding plot lines displaying $h w t$ and $h w$ become more transparent, so that the most transparent green line and gray line show $h w t$ and $h w$ of the 10 dB case, respectively. The comparison of the $x$-axis of Figs. 6(a) and 6(b) indicates that new ordinal partitions occur in the noisy time series. These ten new ordinal partitions are depicted in gray on the $x$-axis of Fig. 6(b) and have negligible $h w t$ and $h w$ amplitude. Figure 6(b) demonstrates that the amplitude of the entropies changes with the inclusion of noise. Despite this, the rank ordering of the ordinal symbols appears to be robust. Specifically, the ranking of $h w$ is nearly the same up to 20 dB and the ranking of $h w t$ is nearly the same up to 40 dB.

Figure 7 shows the entry points for each section more clearly, with larger and non-transparent points both on the time series [Fig. 7(a)] and the embedded attractor [Fig. 7(b)]. These points are the ones that are used to calculate the weight in $h w t$⁠. Each set of dots in each ordinal partition has created an FRM, which is shown in Fig. 7(c). Comparing Figs. 7(c) and 1(c) indicates that the FRM of the descending ordinal partition [4, 3, 2, 1] (light green) gives the same FRM as the local maxima. The FRMs are divided into two wings of the attractor by the black diagonal dashed line in the middle of Fig. 7(c). All the FRMs obtained from other ordinal partitions in the first and second entropy levels [the colors of the first 10 rows of Fig. 5(a)] have a section above and a section below this diagonal and represent the main features of the attractor. While the FRMs of the third entropy level [the colors of the last four rows of Fig. 5(a)] only have some random points on one side of the diagonal, they do not represent any particular information about the topological structure of the attractor.

As previously stated, the attractor of any dynamical system can be divided into different ordinal partitions based on any desired number of points in the ordinal window. As a result, the FRMs derived from these ordinal partitions can represent the attractor’s dynamical features. The difference between the FRMs of different ordinal partitions grows larger as the attractor’s dynamics become more complicated. As shown in the  Appendix, the difference between the FRMs of different ordinal partitions for the Mackay–Glass system is greater than the difference between the FRMs of different ordinal partitions for Lorenz and Rössler as it has much more complicated dynamics.

To investigate the effect of noise on the FRMs, a similar method was applied to the noisy time series. As shown in Fig. 7(c), the FRM of the ordinal partition [4 3 2 1] is very similar to the FRM of local maxima. As a result, to investigate the robustness of the presented method to noise, we compare these sets of FRMs. Figure 8 depicts the FRMs of local maxima (top panel) and the FRMs of the ordinal partition [4 3 2 1] (bottom panel) for noise levels of 100 dB [Figs. 8(a) and 8(e)], 70 dB [Figs. 8(b) and 8(f)], 50 dB [Figs. 8(c) and 8(g)], and 30 dB [Figs. 8(d) and 8(h)]. The comparison of the FRMs of local maxima with 70 dB [Fig. 8(b)] noise level by 100 dB [Fig. 8(a)] which demonstrate almost noise-free signal reveals that three short lines have been added to the bottom right and also top and bottom left of the FRM. However, only two dots have been added to the FRM extracted from the descending ordinal partition in 70 dB [Fig. 8(f)]. It can be seen that the FRMs of descending ordinal partition in 50 dB [Fig. 8(g)] and 30 dB [Fig. 8(h)] are very similar to the FRMs of local maxima in 70 dB [Fig. 8(b)] and 50 dB [Fig. 8(c)], respectively, implying that the FRMs of descending ordinal partition follow the FRMs of local maxima by adding 20 dB more noise. This robustness of FRMs of high-rank ordinal partitions to noise has also been observed in other systems such as Rössler.

We presented a new generic method for analyzing dynamical systems based on their output time series that can be applied to a wide range of dynamical systems regardless of complexity. This method extracts Poincaré sections by using ordinal partitions without the use of an embedded attractor or generating a high-dimensional complex network. The novel types of weighted permutation entropies presented here aid in ranking the effectiveness of these Poincaré sections. Since ordinal partitions with high entropy values have a larger sector on the attractor, they can be considered suitable Poincaré sections for analyzing the system’s dynamical behaviors. The combination of FRMs from all of these sections can provide a group of system information that aids in the analysis. A group of FRMs can either allow one to choose their favorite or provide a more comprehensive picture of the dynamics.

We investigate the application of this method to various dynamical systems and, in all cases, are able to represent several one-dimensional maps of the system accurately, each capturing the system’s main dynamical characteristic. In particular, our method generates various types of FRMs that show detailed aspects of the system’s dynamics in more complicated dynamical systems (such as Mackey–Glass). Embedding demonstrates why all of the high entropy ones are good and capture the system accurately. Some of these good FRMs are comparable to the well-known Poincaré sections (such as $x ˙=0$⁠) with the advantage that they are much more noise resistant.

While the method we have proposed here is directed at the classic construction of FRMs—successive intersections of Poincaré sections—the method itself could be applied without modification to any time series, making our approach superior to traditional methods such as finding maxima. Without the underlying continuity of a flow, we would only be able to look at successive first inclusion within a particular ordinal class. This would still represent a projection to a lower dimensional system, but the classical connection between the flow and the first return map has no obvious immediate analog. This may also yield a valuable technique for a discontinuous time series analysis, but the analytic underpinnings are less obvious (to us).

We wish to acknowledge the Australian Research Council funding for DP200102961. S.D.A. also acknowledges funding from the Forrest Foundation.

The authors have no conflicts to disclose.

Zahra Shahriari: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal). Shannon D. Algar: Conceptualization (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). David M. Walker: Conceptualization (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Michael Small: Conceptualization (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article.

This section provides examples of our method’s application to other dynamical systems.

Figure 9 depicts the results, equivalent to Figs. 6 and 7, for the case when $m=10$⁠.