A Look into Chaos Detection through Topological Data Analysis

TL;DR

This paper introduces a novel chaos detection method using topological data analysis, capable of analyzing noisy time series without prior models, and compares favorably with existing tests on classical chaotic systems.

Contribution

The paper presents a new chaos detection approach based on topological data analysis that works reliably on noisy, model-free time series and distinguishes chaotic from periodic dynamics.

Findings

Effective on noisy data and model-free scenarios

Comparable to the 0--1 correlation test on clean data

Fails to converge on trajectories with partially predictable chaos

Abstract

Traditionally, computation of Lyapunov exponents has been the marque method for identifying chaos in a time series. Recently, new methods have emerged for systems with both known and unknown models to produce a definitive 0--1 diagnostic. However, there still lacks a method which can reliably perform an evaluation for noisy time series with no known model. In this paper, we present a new chaos detection method which utilizes tools from topological data analysis. Bi-variate density estimates of the randomly projected time series in the $p$-$q$ plane described in Gottwald and Melbourne's approach for 0--1 detection are used to generate a gray-scale image. We show that simple statistical summaries of the 0D sub-level set persistence of the images can elucidate whether or not the underlying time series is chaotic. Case studies on the Lorenz and Rossler attractors as well as the Logistic Map are used to validate this claim. We demonstrate that our test is comparable to the 0--1 correlation test for clean time series and that it is able to distinguish between periodic and chaotic dynamics even at high noise-levels. However, we show that neither our persistence based test nor the 0--1 test converge for trajectories with partially predicable chaos, i.e. trajectories with a cross-distance scaling exponent of zero and a non-zero cross correlation.

Figures

Equations

$$\dot{x}=\sigma (y-x),\dot{y}=x(\rho -z)-y,\dot{z}=xy-\beta z.$$

$$C_{12}=\frac{1}{s^2}⟨\left[{\mathbf{x}}_1(t)-E({\text{x}}_1)\right]\cdot \left[({\mathbf{x}}_2(t)-E({\text{x}}_2)\right]⟩,$$

$$E(\mathbf{x})=\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{j=1}^N\mathbf{x}(j),$$

$$s^2=\underset{T-\infty }{\lim }{\int }_T^{2T}{\left[\mathbf{x}(t)-E(\text{x})\right]}^2\mathrm{d}t.$$

$$d_{12}(t>T_{\lambda })\propto {\delta }^{\nu },d_{12}(t)=⟨\left|{\mathbf{x}}_1(t)-{\mathbf{x}}_2(t)\right|⟩$$

$$p(n)=\sum _{j=1}^n\varphi (j)\cos jc,\mathrm{and}q(n)=\sum _{j=1}^n\varphi (j)\sin jc,$$

$$M_c(n)=\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{j=1}^N{\left[p_c(j+n)-p_c(j)\right]}^2+{\left[q_c(j+n)-q_c(j)\right]}^2.$$

$$V_{\mathrm{osc}}(c,n)={\left(E\varphi \right)}^2\frac{1-\cos \left(nc\right)}{1-\cos \left(c\right)}.$$

$$D_c(n)=M_c(n)-V_{\mathrm{osc}}.$$

$$\mathrm{cov}\left(x,y\right)=\frac{1}{q}\sum _{j=1}^q\left(x\left(j\right)-\bar{x}\right)\left(y\left(j\right)-\bar{y}\right),\mathrm{and}\mathrm{var}\left(x\right)=\mathrm{cov}\left(x,x\right),$$

$$K_c=\mathrm{corr}\left(\xi ,\Delta \right)=\frac{\mathrm{cov}\left(\xi ,\delta \right)}{\sqrt{\mathrm{var}\left(\xi \right)\mathrm{var}\left(\Delta \right)}}\in [-1,1],$$

$$K_c=\underset{n\to \infty }{\lim }\frac{\log {\tilde{D}}_c(n)}{\log n},\widetilde{D_c}(n)=D_c(n)-\underset{n=1,\dots ,n_{\mathrm{cut}}}{\min }{D}_c(n).$$

$$K_c=\underset{n\to \infty }{\lim }\frac{\log M_c(n)}{\log n}.$$

$$L_{\lambda }=\left\{x:f\left(x\right)\le \lambda \right\}=f^{-1}\left(\left[-\infty ,\lambda \right]\right).$$

$$\hat{f}(x,\mathrm{H})=\frac{1}{n}\sum _{i=1}^n{\text{K}}_{\mathrm{H}}\left(x-x_i\right),$$

$$G_h=\frac{1}{2\pi h^2}{e}^{-\frac{x^2+y^2}{2h^2}}$$

$$PS_{1}=\Bigg{\langle}\sum_{j=1}^{N_{i}}\left\{\frac{\ \sqrt{(d_{j}^{2}+b_{j}^{2})}}{N_{i}},(b_{j},d_{j})\in{\mathsf{D}}_{i}\right\}\Bigg{\rangle},$$

$$\dot{x}=-y-z,\dot{y}=x+ay,\dot{z}=\gamma +z(x-c),$$

$$x_{n+1}=\mu x_n(1-x_n).$$

Full text

\unmarkedfntext

Address correspondence to khasawn3@egr.msu.edu

\unmarkedfntext

Code and data available at Mendeley repository Chaos Detection with Persistent Homology

\unmarkedfntext

Accepted February 26, 2020 in Physica D: Nonlinear Phenomenon, DOI: 10.1016/j.physd.2020.132446

**A Look into Chaos Detection through Topological Data Analysis

**

Joshua R. Tempelman and Firas A. Khasawneh

(*Department of Mechanical Engineering, Michigan State University, 428 S. Shaw Ln. East Lansing, MI 48824

*)

Abstract. Traditionally, computation of Lyapunov exponents has been the marque method for identifying chaos in a time series. Recently, new methods have emerged for systems with both known and unknown models to produce a definitive 0–1 diagnostic. However, there still lacks a method which can reliably perform an evaluation for noisy time series with no known model. In this paper, we present a new chaos detection method which utilizes tools from topological data analysis. Bi-variate density estimates of the randomly projected time series in the $p$-$q$ plane described in Gottwald and Melbourne’s approach for 0–1 detection are used to generate a gray-scale image. We show that simple statistical summaries of the 0D sub-level set persistence of the images can elucidate whether or not the underlying time series is chaotic. Case studies on the Lorenz and Rossler attractors as well as the Logistic Map are used to validate this claim. We demonstrate that our test is comparable to the 0–1 correlation test for clean time series and that it is able to distinguish between periodic and chaotic dynamics even at high noise-levels. However, we show that neither our persistence based test nor the 0–1 test converge for trajectories with partially predicable chaos, i.e. trajectories with a cross-distance scaling exponent of zero and a non-zero cross correlation.

Key words. Topological data analysis, chaos detection, persistent homology, time series

1 Introduction

Chaos is an obscure phenomenon which appears in many physical and theoretical systems. Chaotic systems display a hyper-sensitivity to initial conditions, and their subsequent dynamics are difficult to anticipate as they are seemingly random and sporadic. Differentiating a chaotic time series from a periodic one is not a trivial task. The traditional recipe for doing so utilizes a tangent space method where, for example, the sign of the maximal Lyuponov exponent is used to classify the systems behavior. This method has proven reliable, although it is not straightforward to implement for systems with no known model, even with the procedures given by Bennetin [1] and Wolf et al. [2]. Additionally, the tangent space methods perform poorly when dealing with noisy finite time series, and this is a major limitation when analyzing data from a high-dimensional and noisy system.

There have been recent efforts to detect chaos in time series such as the binary test introduced by Wernecke et al. [3] which can also distinguish strong chaos from partially predictable chaos. Like the Lyuponov exponent procedure, this method is reliant on the time–evolution of initially close trajectories. Thus, a known model is required to deliver a diagnosis. Additionally, there are no results given in [3] on the effects of noise on the scoring method.

A popular emergent tool for chaos detection in a time series is the 0–1 test which was introduced by Gottwald and Melbourne in 2004 [4]. Unlike tangent space methods, knowledge of the full state–space is not required to perform the 0–1 test.

Since a state-space reconstruction is not required, the 0–1 method is favorable for data streams for which not all state variables are observable; this is the case for many engineered and natural systems.

The 0–1 test has been presented in two forms: the regression test and the correlation test (see Section 1.2 and Refs. [4, 5, 6, 7]). The correlation test has been shown to perform exceptionally well for noise-free time series [6]. However, when noise contaminates the data, a modified version of the original regression test desensitizes the scores to noise [5, 6]. Alternatively, the correlation method handles noise exceptionally well if “gauge” values are given for a time series known to be periodic [8]. However, the required gauge score is not always obtainable in real world systems.

In this paper, we present an alternative approach for combating the sensitivity of the 0–1 test to noise. Our approach combines elements from the 0–1 recipe with a tool from Topological Data Analysis (TDA) specifically persistent homology. We leverage the stability of persistent homology in the presence of noise to deliver a method for identifying the shift-points between chaotic and regular dynamics in moderately noisy time series.

The paper is organized as follows. In Sections 1.1 and 1.2, we introduce the dynamic regimes we are studying and review the 0–1 test, respectively. In Section 2, we describe persistent homology and how this tool can be used in conjunction with projections of the time series data in order to identify the shift point between chaotic and periodic dynamics. Section 3 investigates the robustness of our approach to noise, and the summarizing conclusion is given in Section 4.

1.1 Definitions of the Studied Dynamic Behavior

The exposition of our method uses simulations of the Lorenz model. The equations of motion for this dynamical system are defined as [9]

[TABLE]

where nominal definitions of the Lorenz constants for this study are $\beta=8/3$, $\sigma=10$, and $\rho\in[180.3,\ 181.3]$ is left as a bifurcation parameter. The manipulation of $\rho$ can lead to three distinct dynamic regimes as shown in Fig. 1. In [3], it is shown that the Lorenz system displays periodic, partially predictable chaos (PPC), and chaotic behavior over the range $\rho\in[180.3,\ 181.3]$. These regimes are formally defined in [3] as strong chaos, PPC, and laminar flow. We implement the 0–1 tests of [5, 6] as well as our proposed persistence–based method on time series which were shown in [3] to be PPC. The PPC trajectories are characterized by having a non-zero cross correlation for time-scales beyond the Lyapunov prediction time and a cross-distance scaling exponent $\nu$ of zero. The cross-correlation is defined as

[TABLE]

where $E(\textbf{x})$ is the expected value given by

[TABLE]

$\langle\hskip 3.0pt\rangle$ denotes the ensemble mean, and the variance $s^{2}$ is

[TABLE]

The cross-distance scaling between trajectories $\textbf{x}_{1}$ and $\textbf{x}_{2}$ is defined with

[TABLE]

where $T_{\lambda}$ is the Lyapunov prediction time $T_{\lambda}=\ln\left(|\textbf{x}_{1}-\textbf{x}_{2}|/\delta\right)/\lambda_{m}$, $\delta$ is the distance of initially close trajectories, and $\lambda_{m}$ is the maximal Lyapunov exponent. The reader is directed to [3] for a detailed explanation and demonstration of these definitions.

In this paper, the regimes will be referred to as chaotic, PPC, and periodic dynamics, respectively. The ensemble mean of $10^{4}$ independent trials used in [3] is very effective in identifying these respective regimes with a binary diagnostic (Fig. 1). However, these evaluations require a known set of differential equations to generate neighboring initial conditions. Thus, the approach in [3] requires the full knowledge of the state–space and cannot be applied to a single observation of a high–dimensional system.

1.2 The Traditional 0–1 Test

A 0–1 diagnostic for identifying chaos in time series with known or unknown models, i.e. using recorded data, is given in [4]. The method has been presented in several different forms since its original conception [4, 5, 6, 7]. In this paper, we use the modified regression method described in [5] as well as the correlation method as described in [7]. The original test identifies some observable characteristic $\phi(j)$ of a dynamic system and uses it to construct planar $(p,q)$ data sets according to

[TABLE]

where $c$ is a random variable drawn from a uniform distribution in $(0,\pi)$ and $n=\left(0,1,\dots,{N}\right)$. Using a single data set of length ${N}$, the number of $c$ values selected creates $N_{c}$ data sets in the $(p,q)$ plane. In this paper, $N_{c}=100$ and $N=5000$ unless otherwise stated. The $p$-$q$ data is then used to compute the mean-square displacement

[TABLE]

Note that in Eq. (6), $n=\left(1,2\dots,n_{\rm cut}\right)$ where $n_{\rm cut}=N/10$. The test for chaos is based on the asymptotic growth of $M_{c}(n)$ with respect to $n$, which is a bounded function of time for periodic dynamics and grows linearly with time for chaotic dynamics [6]. It is shown in [6] that the modified mean-square displacement $D_{c}(n)$ possesses better convergence characteristics. To construct $D_{c}(n)$, begin by defining the oscillatory term

[TABLE]

Now, subtract the $V_{\rm osc}$ term from the mean-square displacement to define the modified mean-square displacement

[TABLE]

The test can now be performed using either the correlation test or regression test. Using the standard definitions for covariance and variance respectively

[TABLE]

the correlation $K_{c}$ value is found using the correlation function

[TABLE]

where $\xi=\left(1,2,\dots,n_{\rm cut}\right)$ and $\Delta=\left(D_{c}(1),D_{c}(2),\dots D_{c}\left(n_{\rm cut}\right)\right)$. The 0–1 correlation $K$ value of the time series is then defined as the median value of the $K_{c}$ values with a median near zero noting a periodic time series and a median near one noting a chaotic time series.

Alternatively, the 0–1 regression test may be applied according to

[TABLE]

In [5], a $K_{c}$ value is computed for $M_{c}(n)$ rather than $\tilde{D}_{c}(n)$, and this will be referred to as the modified 0–1 regression test

[TABLE]

For finite data, the $K_{c}$ score is computed as the slope of the line which fits the $\log(M_{c}(n))$ versus $\log(n)$ with the least absolute deviation [6]. The modified 0–1 regression test will be used instead of the standard regression test in this paper since it has been shown to perform better in the presence of noise [5, 6]; noise robustness will play major role in later sections of this paper (section 3).

The 0–1 test typically works better than tangent space methods for detecting chaos. However, some abating factors must be addressed. The test will always return a periodic diagnosis when the data is over-sampled for the case of time-continuous systems [10]. Additionally, false positives are found when periodic time series are contaminated with substantial noise [5]. The issue of oversampling can be resolved by properly sub-sampling the time series according to [6, 10, 11, 12], the effects of noise are more difficult to overcome. A method has been developed which allows the 0–1 test to handle noise very effectively if a gauge-value (i.e. time series with known periodic parameters) is provided [8]. Furthermore, the effect of noise on the 0–1 test has been leveraged to predict the noise level within a periodic time series [7].

If noise is of concern and using a gauge value is not an option, the modified regression score should be used [5, 8]. However, whereas the 0-1 test as described above has a mathematically rigorous underpinning [13], this is not the case for the more noise-robust modification [5].

An interesting intermediate result of the 0–1 method is the projection of the data onto the $p$-$q$ plane. If the dynamics of the time series are periodic, then the resulting projection will have a compact and typically annular topology. In contrast, Fig. 2 shows that chaotic dynamics lead to an irregular scattering of points in the $p$-$q$ space. Therefore, the topology of the $p$-$q$ projections can give insight into the underlying dynamics, as we show in Section 2.

2 A Topological Approach

The motivation behind a TDA approach is that data can have shape, and that shape has a meaning. TDA provides tools which lead to a computational summary of the shape of data.

We showed in Section 1.2 that the 0–1 test produces data projections with meaningful topology. Thus, these projections can be quantified in a topological sense and studied in a new framework.

The ability to quantify the topology of a data structure has led to many breakthroughs in various fields. Specifically, TDA has already made meaningful contributions in time series analysis [14, 15, 16], economics [17], machining dynamics [18, 19, 20, 21], biochemistry [22], and plant morphology [23]. In this section, we extend the application of this powerful computational tool to chaos detection.

Recently, Mittal and Gupta have explored persistence homology as a tool for detecting early bifurcations and chaos in dynamic systems [24]. In [24], persistent homology is used to identify the presence of holes in the phase portraits of deterministic dynamical systems, thus an indicator for bifurcations. The method has shown promising results in [24] even with a signal to noise ratio of 30 dB. However, wavelet filtering is required to pre–process noisy data. Additionally, we note that their method requires one to have access to the whole phase-space or to perform a phase-space reconstruction. The approach we describe here requires neither data pre–processing nor phase-space reconstruction.

We use the Lorenz system as a test model for elucidating our TDA-based approach for chaos detection. Time series which are generated for different values of the bifurcation parameter $\rho$ in the Lorenz model are used to exhibit chaotic, PPC, and periodic dynamics. The Lorenz time series were produced using the ODE45 function in MATLAB® with a relative tolerance of 1e-6, a time step of $\Delta{t}=0.001$, and a transient cut-off of 100 seconds.

Consistent with the 0–1 test, the time series data must first be sub-sampled, if needed. This was done using the frequency approach presented by Melosik and Marszalek [10, 11]. This sub-sampling method results in a time series sampled such that $2f_{\rm max}<f_{s}<4f_{\rm max}$ where $f_{s}$ is the sub–sampling frequency and $f_{\rm max}$ is the maximum significant frequency of the data in the frequency domain. In the 0–1 test, failure to sub-sample may result in false-negatives. Similarly, we found that oversampling affects the topological characteristics of the $p$–$q$ projections thus leading to incorrect conclusions when studying the topology of the $p$–$q$ space. As a general rule-of-thumb, it is suggested that the sampling rate chosen in the time series sub–sampling is set to three times the maximum significant frequency of the data.

The subsequent sections give a detailed development of the TDA approach; the overview is as follows. Multiple projections of $p$-$q$ data are generated for a given time series. The bi-variate Kernel Density Estimates (KDEs) of these projections produce a gray scale image (Fig. 4). Sub-level set persistence is then used to provide a computational summary of the topological structure of the KDE. We give a brief explanation of persistent homology in Section 2.1 and the procedure for applying TDA to the $p$-$q$ data in Section 2.2. The results are explained in Section 2.3 while case studies on the Rossler system and the Logistic map are given in Sections  2.3.1 and 2.3.2 respectively. The implementation of this method is described with Algorithm 1 which we implement in MATLAB® and is published online [25].

2.1 Sub-level Set Persistence

Only a brief explanation of persistent homology, commonly known as persistence, will be provided in this paper. For a more comprehensive explanation, the reader is directed to [26, 27, 28, 29, 30, 31].

Persistence can be applied to 3-dimensional data such as contour plots or even gray-scaled images. Each point in topological space $\mathbb{X}$ is assigned a real number, which for the purpose of this explanation, may be thought of as a height coordinate. This resulting function $f:\mathbb{X}\to\mathbb{R}$ can have local maxima and minima. Persistent homology is interested in changes in topology at various sub-level sets. For this function $f$, the $\lambda$ sub-level set is

[TABLE]

For any pair of levels $\lambda_{1}>\lambda_{2}$ the condition $L_{\lambda_{1}}\supseteq L_{\lambda_{2}}$ must be satisfied and the collection of sets $\{L_{\lambda}\}_{\lambda\in{\mathbb{R}}}$ forms a filtration where the level is the index set [32]. Consider Fig. 3 which shows a toy example for computing 0D sub-level persistence for a function $f$. One way to visualize the process is to imagine a rising water level beginning at the base of $\mathbb{X}$ and account for each time a new minimum of the surface is touched and new maximum is reached (see Fig. 3). As the level $\lambda$ is raised, new topological features are captured through generators of its homology group. Two dimensions of homology groups are considered. The 0D groups capture connected components (see Fig. 3) while the 1D groups captures regions forming circular structures (see Fig. 5).

As new levels are introduced, new generators of homology groups are formed and existing generators merge. The level at which these generators and groupings occur are the birth and death time, respectively, of the topological feature. Therefore, the three characteristics for each generator are the homology dimension, birth time, and death time. The collection of these points is drawn as what is referred to as a persistence diagram. For $|D|$ generators, the persistence diagram is ${\mathsf{D}}=\left\{(r_{j},b_{j},d_{j}):j=1,\dots,|D|\right\}$ where $r_{j}$, $b_{j}$, and $d_{j}$ are, respectively, homology order, birth time, and death time of the $j$th generator [32].

2.2 1D Persistence of Kernel Density Estimates

In this paper, we utilize the 0D sub-level set persistence to extract the underlying characteristics of the probability densities of the $p$-$q$ projections. The KDEs correspond to the $p$-$q$ data according to [33]

[TABLE]

where ${x}=(p,q)$, $n$ is the number of elements in $\{x_{1},x_{2},\dots,x_{n}\}$, and ${\textbf{K}_{\mathrm{H}}}$ is the kernel function. The resulting KDEs are converted to a gray scale image. The gray scale image is normalized to have a pixel value in [0, 1].

It should be noted that increasing the grid resolution will drastically lengthen the runtime of the program. To reduce computation time of the corresponding sublevel set persistence, a KDE grid size of $2^{5}\times 2^{5}$ is used. This is a relatively coarse grid resolution and results in KDEs which are jagged when compared to a KDE of a higher resolution. This could lead to misleading computational summaries of the KDE topology, thus a 2D Gaussian smoothing function

[TABLE]

is used. We found $h=1.3$ to be a good choice for the smoothing kernel width at a $2^{5}\times 2^{5}$ grid resolution (the effect of this kernel width on the results of the persistence test is discussed in Section 3.1). The DIPHA111https://github.com/DIPHA/dipha library was used to generate all the needed 0D persistence diagrams.

2.3 Results and Discussion

The central motivation for this study is to yield a reliable interpretation of the underlying dynamics of time series data using persistence. To better understand how this may be achieved, refer to Fig. 6. Here, multiple persistence diagrams are overlaid for ranges of the bifurcation parameter $\rho$ of the Lorenz model. The persistence points with a birth time near zero are omitted from Fig. 6 since they are considered noise in this application. The persistence points shown in Fig. 6 show that when persistence diagrams are generated from time series associated with a periodic bifurcation parameter, the persistence points accumulate in the upper-right corner of the persistence diagram.

These produce birth/death combinations which are tightly clustered and bounded to regions on or near the diagonal. There is an apparent correlation here between the birth and death times of the persistence points. Furthermore, these points all appear to be occurring at higher birth and death times than what is seen for their chaotic counterparts. The persistence data generated form PPC trials are visually indistinguishable from periodic trials. On the other hand, chaos based persistence diagrams are visually perceptible from periodic when many trials are overlaid as shown in Fig. 6. Thus, the hypothesis to test is whether or not this behavior may be reliably used to trace back the dynamics of the original time series.

Scoring: Infinitely many $p$-$q$ projections may be generated for a given time series since $c$ may be indefinitely drawn from $(0,\pi)$. Thus, we use an ensemble of $p$-$q$ projections and use a statistical summary of birth and death times to deduce the original dynamics of the time series. In this section, the mean distance of persistence points from the diagram’s origin are used to discern the underlying dynamics of the time series.

The threshold for chaotic versus periodic or PPC behavior in persistence diagrams can be defined as a result of the location of the points in the diagram. A basic statistical evaluation of the trials used in Fig. 6 is used to elucidate how such a threshold is defined. Looking at the average birth and death times, Fig. 7 shows that the average time of birth and death for periodic data is substantially higher than the chaotic counterpart (approximately $0.90$–$1.0$ for periodic and $0.40$–$0.65$ for chaotic). The statistical summary outlined in Fig. 7 ignores the birth/death combination on or near zero birth time vertical. These persistence points do not reveal any relevant information, and can be thought of as irrelevant noise in the KDE plots.

The statistical summary was conducted by taking the mean birth and death time of each diagram for 200 realizations of $p$-$q$ projections at each bifurcation parameter and then computing the ensemble mean of these values. The average birth and death times are substantially higher for periodic and PPC data; therefore, we score the data using the average distance to the origin according to

[TABLE]

where $PS_{1}\in[0,\sqrt{2}]$, ${\mathsf{D}}_{i}$ is the persistence diagram generated form the $i$th selection of $c\in(0,\pi)$, and $N_{i}$ is the number of 0D persistence points in ${\mathsf{D}}_{i}$.

The reason for the higher average birth and death times in the persistence diagrams for KDEs generated from periodic time series as compared to chaotic ones is illustrated in Fig. 8. The diffusive nature of the $p$-$q$ trajectories [8] for chaotic data results in KDEs which typically have multiple local maxima and minima randomally appearing in the $p$-$q$ space. These local extrema correlate to the births and deaths of homology classes in the persistence diagram. Alternatively, the bounded projections of periodic time series data will consistently produce a KDE which is bound to a finite region in the $p$-$q$ space with prominent peaks representing frequently visited parts of the $p$-$q$ space. This structure of the KDE leads to the homology classes being constant during the filtration until the top of the structure is reached. At which point, ridges near the peak of the KDE produce many generators near a value of one which die soon after birth.

Thus, it can be concluded that for the Lorenz system, if the $PS_{1}$ score for a set of $p$-$q$ projections is near $\sqrt{2}$, then the underlying dynamics are periodic. The mean distance from the origin $PS_{1}$ was consistently found to be a successful score for classifying chaos for the systems studied in this paper and thus it will be the focus of the results section.

Figure 9 shows a direct comparison between our proposed method (see Fig. 4 and Algorithm 1) and the 0–1 correlation test. The 0–1 test and $PS_{1}$ scores are performed using $N=5000$ (after sub-sampling), and a long time limit ($N=25000$) evaluation of the 0–1 test is included as a reference case [8, 34]. When the dynamics are periodic ($\rho>180.73$), the distance from the origin stabilizes near 1.3. When chaotic dynamics are present, the $PS_{1}$ score stabilizes at approximately 0.8. The scores of 1.3 for periodic and 0.8 for chaotic were found to be consistent when testing with longer time series as well.

The numerous numerical experiments that we preformed seemed to yield $PS_{1}$ scores close to the $[1/\sqrt{2}\approx 0.7$, $\sqrt{5/9}\approx 0.75]$. These two scores coincide with the geometric centers of the diagonal line and the upper left triangle of the persistence diagram, respectively. This suggests that the persistence of the corresponding probability density estimates of the projections possess a global shape which appears similar to the one obtained form the 0D persistence of Gaussian random field excursions [35]. However, computing distributions of persistence diagrams is not a trivial task [36, 37, 38, 35] and is complicated by factors such as the highly nonlinear nature of the resulting metric space, as well as the non-uniqueness of geodesics and means in the persistence diagram space. Finding a closed-form expression for these expectations is a topic for future research.

Figure 9 also shows how the 0–1 test and the $PS_{1}$ scores behave when PPC dynamics are present in the bifurcation parameter range $\rho\in[180.7,181]$. Specifically, in this parameter range both methods return a periodic diagnostic for $N=5000$. However, Fig. 10 shows that the long-time limit ($N=10^{6}$ points) for PPC data deviates from this periodic result for both the 0–1 correlational and regression scores. The figure shows that as the time series is extended, both the 0–1 correlation test and the $PS_{1}$ statistic start to yield non-binary results between zero and one, but not necessarily near zero or one. Although this is found to be consistent in both the modified regression and correlation 0–1 tests, the divergence from a periodic score occurs quicker with the latter. Similarly, a $PS_{1}$ score between 0.8 and 1.3, but not necessarily near 0.8 or 1.3, is returned from the $PS_{1}$ statistic. One may argue that this is simply an issue of finite time. While it appears that the PPC data are tending towards one of the two diagnostic limits for both the 0–1 tests and $PS_{1}$ scores, it does so at too slow of a rate for convergence to be achieved without resorting to excessively taxing computations. In contrast, scores for strongly chaotic or fully periodic regimes remain consistent during this long-time evaluation. This shows that it is difficult to distinguish PPC solely from $p$-$q$ projections.

2.3.1 Case Study: The Rossler System

In order for our method to be applicable to a wide class of signals, it must yield consistent numerical results for periodic and chaotic time series regardless of the source of the signal. To investigate this, we use the Rossler system to gather additional synthetic data. Like the Lorenz equations, the Rossler model is a classic system which has been studied widely [39] and it is described by

[TABLE]

where $a$, $c$ and $\gamma$ are scalars. We only consider fully periodic and fully chaotic signals in this section. The Rossler data is generated through numerical simulation using the MATLAB® ODE45 integration tool with $\Delta{t}=0.01$, a relative tolerance of $10^{-5}$, and a transient cut-off of 1,000 seconds. The parameters selected for this study are $a\in[0.25\hskip 4.0pt0.55]$, $\gamma=2$, and $c=4$ [40]. Figure 11 shows the results of the 0–1 correlation test along with the $PS_{1}$ scores for $N=5000$ as well as the long time reference using $N=25000$ using the 0–1 correlation test. The bifurcation diagram and trajectories in the phase space are also given in Fig. 11a and 11b, respectively.

The same TDA based approach we used for the Lorenz system for chaos detection is applied to the Rossler data. Time series were generated for a span of bifurcation parameters. The minimum significant frequency sub-sampling method described in section 2 (i.e. $2f_{\rm max}<f_{s}<4f_{\rm max}$) was applied to the Rossler simulations which are then projected into the $p$-$q$ space per Eq. (5). The 0D sub-level set persistence then generated the birth/death times for gray-scale renditions of the normalized KDE of these projections which were smoothed per Eq. (15) with $h=1.3$.

The $PS_{1}$ and the correlation 0–1 scores were obtained for the time series generated over the different values of $a$ with $N=5000$. Figure 11 shows that the correlation 0–1 test and the $PS_{1}$ scores perform comparably for $N=5000$, with a more definitive classification $PS_{1}$ scores for $a\in[0.375,0.425]$. In this range of $a$ values, the regression 0–1 test values are rather ambiguous since the $K$ values commonly fall between 0.5 and 0.8 over these selections of $a$.

On the other hand, the $PS_{1}$ and 0–1 scores do not consistently yield a definitive periodic or chaotic diagnostic for $a\in[0.4,0.55]$. While there is a “gray area" in this parameter range, Fig. 11 shows the emergence of elevated plateaus in the $PS_{1}$ scores and 0–1 correlation scores for $a>0.38$. These values coincide with windows of periodicity shown in the bifurcation diagram and in the long–time reference scores of the 0–1 correlation test. Such gray areas are to be expected when dealing with finite time series since the length of the time series may not be sufficient to capture the full dynamics of the attractor [8, 34].

2.3.2 Case Study: The Logistic Map

We further test our method on the Logistic map which is an iterative map governed by [41]

[TABLE]

This nonlinear map was introduced as a simplified population model where $x_{n}$ is a ratio between existing population and maximum population. Here, the bifurcation parameter is $\mu$ with transitions in the dynamics occurring during $\mu\in[3,4)$.

Both the 0–1 correlation test and the persistence scoring method are applied to time series generated from Eq. (18) over the range of $\mu$ values. Like before, the kernel width $h$ in Eq. (15) is set to 1.3 and the persistence points found very close to the origin are ignored. The 0–1 correlation test does a very good job of correctly separating regions of chaos and periodicity through nearly perfect binary scores. The $PS_{1}$ scores do not separate as cleanly as shown in Fig. 12. However, the transition from periodic behavior to chaotic behavior ($\mu\approx 3.55$) as well as the pockets of periodicity in $\mu\in[3.55,4)$ are clearly noticeable in the persistence scores. Additionally, Fig. 12 shows that $PS_{1}\approx 1.4$ indicates a periodic time series while $PS_{1}\approx 0.80$ indicates chaos in the Logistic map. This is consistent with the results in Fig. 9 for the Lorenz system and in Fig. 11 for the Rossler system.

However, it is clear from Figs. 9,  11, and 12 that at times the topological approach does not separate time series data as cleanly as the 0–1 correlation test. Specifically, the transition from chaotic to periodic behavior is at times less obvious in comparison to the original 0–1 test (see Figs. 11 and 12).

2.3.3 Alternative Parameters of the Lorenz Model

In addition to testing across multiple systems, we also explore the Lorenz model of Section 1.1, but this time with the system parameters $\sigma=10$, $\beta=0.1$, and $\rho\in[135,\ 150]$. Again, the time series are sub-sampled per the minimum significant frequency method described in section 2 so that $2f_{\rm max}<f_{s}<4f_{\rm max}$. While the periodic scores for $PS_{1}$ are consistent with the other models (approximately $1.4$), the $PS_{1}$ scores for chaotic time series settle near 1.05. On the other hand, while 0-1 correlation test separates this data in a binary fashion, the 0–1 regression test does not. Figure 13 shows the bifurcations over this span of $\rho$ along with the 0–1 and $PS_{1}$ scores. Again, the extended time 0–1 correlation score ($N=25000$) is used as a reference score.

Note that the weak chaos exhibited by the Lorenz model in this parameter range is not to be confused with PPC. An application of the test given in [3] on the full state-space of the model finds that for all instances of $\rho\in[135,\ 150]$, the cross-correlation $C_{12}$ and scaling exponent $\nu$ either both go to zero for $t>T_{\lambda}$ or both retain a finite nonzero value. Therefore, the model is either fully chaotic or fully periodic in this range by the definitions in [3].

With this we conclude that the $PS_{1}$ scores can be less sensitive to weak chaos than the 0–1 correlation test. Nevertheless, the $PS_{1}$ scores still clearly identify the shift point between periodic and chaotic regimes. This decreased sensitivity of the $PS_{1}$ score is balanced by its increased robustness to noise in periodic time series as we show in Section 3. It is also worth stating that the modified 0–1 regression test (the test prescribed to handle noisy data) performs poorly in this parameter range as well; the $PS_{1}$ scores are in fact more informative here since they yield consistent results for $\rho\in[135,140]$ and $\rho\in[145,150]$ while the $K$ scores of the 0–1 regression test do not.

3 Robustness to Noise

All of our results thus far have addressed simulated, noise free time series. However, all signals obtained in the physical world are inherently contaminated with noise to some degree, and it is important to consider the effect of this noise on new time series analysis methods. Thus, we investigate the robustness of the $PS_{1}$ scores to noise and compare their performance to both the correlation and regression 0–1 methods.

While the correlation 0–1 test does an excellent job in distinguishing chaos from periodicity, the increased sensitivity makes it more susceptible to noise [6]. A brief noise sensitivity demonstration is given in Fig. 14 to exhibit this for simulations of the Lorenz system. This is done by averaging the median $K_{c}$ values as computed by Eq. (10) for time series containing signal-to-noise ratios (SNRs) ranging from 60 to 10 decibels (dB). The average $K$ at each SNR is computed using 100 realizations of these noisy time series. Figure 14 shows that low noise levels (i.e. high SNR) do not affect the resulting K score of the 0–1 correlation test. However, at an SNR of 40 dB, the test begins to yield an ambiguous (non-binary) result. Noise levels above an SNR of 30 dB consistently produce a false-positive indicator for periodic time series. These negative effects on the 0–1 test in the presence of noise were discussed in [5, 34], These references recommend using the 0–1 regression test with noisy time series which was shown to outperform tangent space methods for chaos detection in noisy signals [5, 6].

We start by showing the effect of increasing the level of noise on the structure of the persistence diagrams. Figure 15 shows that these diagrams remain stable as the SNR is varied between 45 and 30 dB. The reason that persistent homology provides this stability is its ability to emphasize the main features of a data set and ignore minor perturbations [28]. In other words, the topological structure of the data does not typically change in the presence of moderate noise. However, that structure starts to break down once noise levels reach the point at which the salient topological features of the KDE are no longer distinguishable from noise.

The robustness of the persistence–based $PS_{1}$ score to noise is further illustrated with Fig. 16. Figure 16 shows the effects of noise on the respective detection methods for the Lorenz model, Rossler model, and Logistic map for SNR levels of 40 and 30 dB. We note that a minimum SNR of 15 dB is often necessary to extract any useful information from a signal. The results obtained from the $PS_{1}$ scores as well as the modified 0–1 test appear stable for 40 dB of noise; however, deviations from the nominal noise-free cases are apparent when the noise levels reach 30 dB. In contrast, the correlation 0–1 test appears to be greatly affected by the noise showing that the trade off of the increased sensitivity to the test is the inability to effectively handle noise.

For the Lorenz model, the last column of Fig. 16 shows that the $PS_{1}$ scores are both stable and binary when noise is added to the signal. While the modified regression method performs well, it does not appear to separate as cleanly as the $PS_{1}$ scores. Similar results are shown for the Rossler model; the correlation method quickly separates from the nominal case when noise levels increase while the modified regression method and the $PS_{1}$ scores remain relatively unaffected. For both the Rossler model and the Logistic maps, the 0–1 correlation test indicates false positives for wide regions of periodic parameters while the persistence scores do not deviate dramatically from what is found in the noise-free data. The $PS_{1}$ scores are comparable to the modified regression 0–1 scores when examining the results of the Logistic map as well. An interesting observation to note is that the correlation method yields a false positive for the periodic parameter $\mu$ near 3.0 even for small noise levels.

3.1 Effect of the Kernel width

The kernel width $h$ of Eq. (15) plays an important role in the final persistence score. The value of $h=1.3$ is selected based on a parameter study of the systems across their bifurcation parameters whereby the value of $h$ is varied between 0.1 and 4. The resulting persistence scores are plotted to reveal which kernel parameter provides the best qualitative separation between periodic and chaotic parameter ranges for both clean and noisy (SNR = 30 dB) time series. Figure 17 shows the effect of the $h$ on the persistence scores. Note that Fig. 17 displays the normalized persistence score $\widetilde{PS}_{1}$, which is normalized to have a maximum of one and minimum of zero.

For both the Lorenz model and Logistics map, there is much better contrast between $PS_{1}$ scores generated from periodic time series versus chaotic time series when $h<2.0$. For the Rossler model, the cleanest separation of chaotic and periodic regimes occurs at $1<h<3.5$. For the Logistics map, the kernel width does not appear to be a factor for the $PS_{1}$ scores ability to appropriately differentiate the classes when no noise is present. However, when the SNR is 30 dB the best range for separation is $1<h<2.5$. Thus, it is concluded that for the systems explored in this paper any value $1.0h\leq 2.0$ will sufficiently distinguish between chaos and periodicity for both clean and noisy time series. Note that if a finer grid resolution is chosen for the KDE construction, this optimal bandwidth range is likely to change.

4 Conclusion

Gottwald and Melbourne’s 0-1 test for chaos uses a one dimensional time series to drive a two-dimensional system in a $p$-$q$ space. The test is based on showing that when the original system’s dynamics are regular, the resulting $p$-$q$ trajectories are typically bounded whereas chaotic dynamics lead to diffusive trajectories. Gottwald and Melbourne used these observations to define two flavors of the test: the 0-1 correlation test, and the 0-1 modified regression test.

The former is a more sensitive test but it incurs larger errors in the presence of noise, while the latter provides better results in the presence of noise. However, mathematical justification has been provided for the 0–1 correlation test [13] while the same cannot be said for the modified 0–1 regression test.

This work builds upon the results of Gottwald and Melbourne and identifies the dynamics of the system that generated the time series by examining the geometry of the density of the $p$-$q$ trajectories using persistent homology. Specifically, A Kernel Density Estimate (KDE) is obtained using the $p$-$q$ projections and the corresponding 0D sub-level set persistence is computed and summarized in persistence diagrams. Several projections are obtained for the same time series to acquire ensemble of persistence diagrams that encode the shape of the corresponding KDEs. Based on observations of 0D persistence on normalized KDEs, we define the persistence score $PS_{1}\in[0,\sqrt{2}]$ which represents the mean distance from the origin for the points in the persistence diagram. $PS_{1}$ assumes values close to $\sqrt{2}$ when the time series is periodic, but its values are below $\sqrt{2}$ for chaotic signals.

The KDEs are smoothed by a Gaussian filter with kernel width $h$. Fig. 17 shows the parameter study which leads to the bandwidth selection of $1.0\leq h\leq 2.0$, and it is found that the results of $PS_{1}$ scores are affected by the kernel bandwidth.

The 0–1 regression and correlation test as well as the $PS_{1}$ scores are applied to the PPC dynamics found in the Lorenz model for $\rho\in[180.7,181,1]$. While Fig. 9 shows that each of these tests report PPC data as periodic in the short-time evaluation ($N=5000$), Fig. 10 shows that each test fails to converge on a binary results for PPC data in the long-time limit ($N=10^{5}$).

For noise-free time series, Figs. 9, 11, and 12 shows that our $PS_{1}$ score successfully delineates between chaotic and periodic behavior albeit not with as clear of a separation as the 0-1 correlation test. In contrast, Fig. 16 shows that when noise is added, the $PS_{1}$ score maintains a consistent performance owing to the stability of persistence diagrams. In this case, comparisons to the more noise robust 0-1 modified regression test show that $PS_{1}$ provides clearer separation especially for high levels of noise, or equivalently, low levels of Signal to Noise Ratio (SNR). If a periodic parameter is known, a gauge value may be obtained and then the test described in [8] can be performed. Otherwise, either the modified regression test or the $PS_{1}$ scores should be evaluated. The $PS_{1}$ scores perform just as well the modified regression test in the presence of noise, and the $PS_{1}$ scores separate data effectively for clean time series.

Acknowledgment

The authors would like to acknowledge Dr. Elizabeth Munch for the helpful insight she provided regarding persistent homology and its applications. This material is based upon work supported by the National Science Foundation under Grant Nos. CMMI-1759823 and DMS-1759824 with PI FAK.