Delay Parameter Selection in Permutation Entropy Using Topological Data Analysis

TL;DR

This paper introduces a topological data analysis method using persistent homology to automatically select the delay parameter in permutation entropy, improving the process's objectivity and reliability across various dynamical systems.

Contribution

It presents a novel TDA-based approach for automatic delay parameter selection in permutation entropy, addressing the challenge of parameter choice in complexity analysis.

Findings

TDA method effectively identifies suitable delay parameters.

Comparison shows TDA approach aligns with expert and optimized parameters.

Method improves objectivity in permutation entropy analysis.

Abstract

Permutation Entropy (PE) is a powerful tool for quantifying the complexity of a signal which includes measuring the regularity of a time series. Additionally, outside of entropy and information theory, permutations have recently been leveraged as a graph representation, which opens the door for graph theory tools and analysis. Despite the successful application of permutations in a variety of scientific domains, permutations requires a judicious choice of the delay parameter $\tau$ and dimension $n$. However, $n$ is typically selected within an accepted range giving optimal results for the majority of systems. Therefore, in this work we focus on choosing the delay parameter, while giving some general guidance on the appropriate selection of $n$ based on a statistical analysis of the permutation distribution. Selecting $\tau$ is often accomplished using trial and error guided by the expertise of domain scientists. However, in this paper, we show how persistent homology, a commonly used tool from Topological Data Analysis (TDA), provides methods for the automatic selection of $\tau$. We evaluate the successful identification of a suitable $\tau$ from our TDA-based approach by comparing our results to both expert suggested parameters from published literature and optimized parameters (if possible) for a wide variety of dynamical systems.

Figures

Tables

Table 1. Table 1: A comparison between the calculated and suggested values for the delay parameter τ 𝜏 \tau . The shaded (red) cells highlight the methods that failed to provide a close match to the suggested delay.

			Delay					Dim.
			1-D Pers. Homol.		Sublevel Set Pers.			R		Exp. Sugg. Parameters			Opt. Param. Range
Cat.	system	State	s	L	t	f	MI	0.5	0.1	$\tau$	$𝐧$	Ref.	$\tau$	$𝐧$
	Gauss.	-	1	1	1	1	3	7	8	1	3-6	[40]	-	-
	Uniform	-	1	1	1	1	3	7	8	-	-	-	-	-
	Rayleigh	-	1	1	1	1	2	7	8	-	-	-	-	-
Noise Models	Expon.	-	1	1	1	1	2	7	8	-	-	-	-	-
		Per.	13	11	11	7	11	5	6
	Lorenz	Cha.	12	13	12	9	12	5	7	10	5-7	[40]	8-17	5-10
		Per.	10	10	10	8	11	5	6
	Rossler	Cha.	12	12	12	10	12	5	6	9	6	[47]	9-15	6-10
		Per.	19	17	16	9	15	5	6
	Bi-direct. Rossler	Cha.	18	16	16	15	17	5	6	15	6-7	[40]	11-22	6-10
		Per.	7	7	6	3	8	5	6
	Mackey Glass	Cha.	7	7	7	4	9	5	7	10	4-8	[50]	6-12	4-8
		Per.	16	17	17	11	19	5	6
	Chua Circuit	Cha.	37	52	17	19	19	5	7	20	5	[44]	16-24	5-10
		Per.	8	8	8	7	9	4	6
	Coupled Ross.-Lor.	Cha.	12	10	8	5	10	5	7	8	3-10	[45]	5-11	4-9
		Per.	16	16	17	11	18	4	5
Cont. Flows	Double Pendul.	Cha.	13	12	10	8	14	6	7	-	-	-	8-20	5-10
	Periodic	-	12	12	13	24	16	4	5	15	4	[47]	-	-
Period. Funct.	Quasi	-	45	46	25	49	26	6	7	-	-	-	-	-
		Per.	1	1	1	1	3	4	5
	Logistic	Cha.	1	1	1	1	16	4	6	1-5	4-7	[40]	1-4	3-6
		Per.	2	2	1	1	3	4	5
Maps	Henon	Cha.	1	1	1	1	16	6	7	1-2	2-16	[40]	1-5	5-8
		Cont.	9	9	22	7	17	5	6
	ECG	Arrh.	13	13	15	6	15	5	6	10-32	3-7	[27]	6-23	5-7
		Cont.	19	18	1	3	6	8	8
Med. Data	EEG	Seiz.	10	4	12	4	10	5	7	1-3	3-7	[40]	2-6	4-7

Equations

$$H(n)=-\sum p({\pi }_i)\log p({\pi }_i),$$

$$H(3)=-\frac{2}{5}\log \frac{2}{5}-\frac{2}{5}\log \frac{2}{5}-\frac{1}{5}\log \frac{1}{5}=1.522bits.$$

$$\begin{array}{rl}{H}_{\max }(n)& =-\sum p({\pi }_i)\log p({\pi }_i)\\ & =-\log \frac{1}{n!}=\log n!.\end{array}$$

$$h_n=-\frac{1}{{\log }_{2}n!}\sum p({\pi }_i){\log }_{2}p({\pi }_i).$$

$$K(G)=\{\sigma \subseteq V\mid uv\in E\text{ for all }u\ne v\in \sigma \}.$$

$${\partial }_n([v_0,\dots ,v_n])=\sum _{i=0}^n(-1{)}^i[v_0,\dots ,{\hat{v}}_i,\dots ,v_n],$$

$$\dots \stackrel{{\partial }_{n+1}}{\to }{\Delta }_n(K)\stackrel{{\partial }_n}{\to }\dots \stackrel{{\partial }_1}{\to }{\Delta }_1(K)\stackrel{{\partial }_0}{\to }0,$$

$$\varnothing \subseteq K(y_1)\subseteq K(y_2)\subseteq \cdots \subseteq K(y_{\ell })=K.$$

$$K(y)=\{\sigma \in K\mid f(\sigma )\le y\},$$

$$f(\sigma )=\underset{\{u,v\}\subseteq \sigma }{\max }d(u,v),$$

$$f({\sigma }_i)=\{\begin{array}{ll}{\chi }_i& \text{ if }{\sigma }_i\in V,\\ \max (u,v)& \text{ if }{\sigma }_i=uv\in E.\end{array}$$

$$K_1\subseteq K_2\subseteq \cdots \subseteq K_N$$

$$H_n(K_1)\to H_n(K_2)\to \cdots \to H_n(K_N).$$

$$I(X;Y)=\sum _{x\in X}\sum _{y\in Y}p(x,y)\log \frac{p(x,y)}{p(x)p(y)},$$

$$S{W}_{m,\tau }f(N)(t)=[f(N)(t),f(N)(t+\tau ),\dots ,f(N)(t+m\tau )],$$

$$f(N)(t)=\frac{1}{|X|}\sum_{k=0}^{N}\Big{(}\sum_{j=0}^{|X|}X(j)e^{-2\pi ijk/T}\Big{)}e^{2\pi ikt/T},$$

$$\ell_{2}(N)=\frac{\Big{(}\sum^{|X|}_{j=0}\big{[}X(j)-f(N)(j)\big{]}^{2}\Big{)}^{1/2}}{\Big{(}\sum^{|X|}_{j=0}X^{2}(j)\Big{)}^{1/2}}.$$

$${\tau }_{\max }=\frac{f_s}{2\min (f_{\mathrm{sig}})},$$

$$\mathcal{L}=\max (\mathrm{pers}({\tilde{\mathcal{D}}}_1)$$

$$s=1-\frac{r_B^2-r_D^2}{3},$$

$$\tau =\frac{f_s}{\alpha f_{\max }},$$

$$z=\frac{x-\mu }{\sigma },$$

$$\mathrm{MAD}=\mathrm{median}(|\mathbf{x}-\tilde{x}|),$$

$$z_m=0.6745\frac{\mathbf{x}-\tilde{x}}{\mathrm{MAD}},$$

$$\frac{N_{\pi }}{n!}\le R,$$

$$\begin{array}{rl}{x}_1(t)& =\frac{t}{10}\\ x_2(t)& =\sin (t)\\ x_3(t)& =\sin (t)+\sin (\pi t)\\ x_4(t)& \sim \mathcal{N}(\mu =0,{\sigma }^2=1),\end{array}$$

$$\Delta h_n(\tau )=h_n^{(Ch.)}(\tau )-h_n^{(Pe.)}(\tau ),$$

$$\mathrm{SN}{\mathrm{R}}_{\mathrm{dB}}=20{\log }_{10}\left(\frac{A_{\mathrm{signal}}}{A_{\mathrm{noise}}}\right),$$

$$A_{\mathrm{noise}}=A_{\mathrm{signal}}{10}^{-\frac{\mathrm{SN}{\mathrm{R}}_{\mathrm{dB}}}{20}}.$$

$$A_{\mathrm{signal}}=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left[x(t_i)-\bar{x}\right]}^2},$$

Full text

1 Introduction

Shannon entropy, which was introduced in 1948 [43], is a summary statistic measuring the regularity of a dataset. Since then, several new forms of entropy have been popularized for time series analysis. Some examples include approximate entropy [37], sample entropy [39], and permutation entropy (PE) [6]. While all of these methods measure the regularity of a sequence, PE does this through the motifs or permutations found within the signal. This allows relating PE to be related to predictability [16, 34], which is useful for detecting dynamic state changes. Similar to Shannon Entropy, PE [6] is quantified as the summation of the probabilities of a data type (see Eq. (1)), where the data types for PE are permutations (see Fig. 1), which we represent as $\pi$. Permutations have recently been used in other applications such as ordinal partition networks [28, 32] and the conditional entropy of these networks [42]. The permutation parameters $n$ and $\tau$ represent the permutation size and spacing, respectively. More specifically, $\tau$ is the embedding delay lag applied to the series and $n$ is a natural number that describes the dimension of the permutation. In this study we focus on selecting $\tau$ using methods based in Topological Data Analysis (TDA) since the dimension is typically chosen in the range $3<n\leq 7$ for most applications [40]. However, we still provide a novel and simple guidance on the automatic selection of $n$ based on a statistical analysis of the permutations.

Currently, the most common method for selecting PE parameters is to adopt the values suggested by domain scientists. For example, Li et al. [26] suggest using $\tau=3$ and $n=3$ for electroencephalographic (EEG) data, Zhang and Liu [48] suggest $\tau=3$ or $5$ and $n\in[3,5]$ for logistic maps, and Frank et al. [14] suggest $\tau=2$ or $3$ and $n\in[3,7]$ for heart rate applications. One main disadvantage of using suggested parameter values for an application is the high dependence of PE on the sampling frequency. As an example, Popov et al. [38] showed the importance of considering the sampling frequency when selecting $\tau$ for an EEG signal. Another limitation is the need for application expertise in order to determine the needed parameters. This can hinder using PE in new applications that have not been sufficiently explored. Consequently, there is a need for an automatic, application-independent parameter selection algorithm for PE.

Several methods have been developed for estimating $\tau$ for phase space reconstruction via Takens embedding [46]. Some of which include mutual information [15], autocorrelation [18], and phase space methods [8]. There has also been recent work in determining if these methods are suitable for the delay parameter selection for permutations [31]. Outside of this, a general framework for selecting both $n$ and $\tau$ was introduced in [40]. In this manuscript our goal is to determine other, TDA-based methods for selecting $\tau$ for PE and to draw a connection between permutations and state space reconstruction. Specifically, in Section 3.2, we relate permutations to state space reconstruction to provide a justification for using the lag parameter from the latter to select $\tau$ for the former. We then present three novel TDA-based tools for finding $\tau$. Our first approach utilizes Sliding Windows and $1$-D Persistence Scoring (SW1PerS) [35] to measure how periodic or significant the circular shape is for the embedded point cloud. Using SW1PerS, we search for the delay that maximizes a periodicity score, which is suggested as an embedding delay for PE due to the connection between phase space reconstruction and permutation entropy as described in Section 3.2. For the second approach we compute the [math]-D sublevel set persistence in both time and frequency domains to obtain approximations of the maximum significant frequency. We then utilize Nyquist’s sampling theorem to find an appropriate $\tau$ value.

To determine the viability of our methods, PE parameters are generated and compared to expert suggested values, optimal parameters based on maximizing the difference between permutation entropy for periodic and chaotic signals, and the delay corresponding to the first minima of mutual information. The method of mutual information was chosen as a basis for comparison based on its accuracy in selecting $\tau$ for PE as demonstrated in [31].

The paper is organized as follows. In Section 1.1 we describe PE and show its computation using a simple example. Section 2 provides an overview of the tools that we use from TDA. In Section 3, we review the method of mutual information for embedding delay parameter selection (Section 3.1), the SW1PerS approach for finding $\tau$ (Section 3.3), and the [math]-D sublevel set persistence methods for selecting $\tau$ (Section 3.4). In Section 4 we introduce our novel approach for selecting the permutation dimension $n$ using the statistics of permutations and highlight the comparison to standard Takens’ embedding techniques which, from our analysis, are not appropriate for permutation dimension selection. The results of each method for a variety of systems are then presented in Section 5, and the concluding remarks are presented in Section 6.

1.1 Permutation Entropy Example

Permutation entropy $H(n)$ for permutation dimension $n$ is calculated according to [6] as

[TABLE]

where $p(\pi_{i})$ is the probability of a permutation $\pi_{i}$, and $H(n)$ has units of bits when the logarithm is of base $2$. The permutation entropy parameters $\tau$ and $n$ are used when selecting the permutation size: $\tau$ is the number of time steps between two consecutive points in a uniformly sub-sampled time series, and $n$ is the permutation length or motif dimension. Using a real-valued dataset $X$ and a measurement of the set $x_{i}\in X$, we can define the vector $v_{i}=[x_{i},x_{i+\tau},x_{i+2\tau},\ldots,x_{i+(n-1)\tau}]$, which has the permutation $\pi_{i}$. To better understand the possible permutations, consider an example with third degree ($n=3$) permutations. This results in six possible motifs as shown in Fig. 1.

Next, to further demonstrate PE with an example, consider the sequence $X=[4,7,9,10,6,11,3]$ with third order permutations $n=3$ and time delay $\tau=1$. The sequence can be broken down into the following permutations: two $(0,1,2)$ permutations, one $(1,0,2)$ permutation, and two $(1,2,0)$ permutations for a total of 5 permutations. Applying Eq. (1) yields

[TABLE]

The permutation distribution can be visually understood by illustrating the probabilities of each permutation as separate bins. To accomplish this, Fig. 2 was created by taking the same series $X$ (Fig. 2a) and placing the abundance of each permutation into its respective bin (Fig. 2b).

PE is at a maximum when all $n!$ possible permutations are evenly distributed or, equivalently, when the permutations are equiprobable with $p(\pi_{i})=p(\pi_{1}),p(\pi_{2}),\ldots,p(\pi_{n!})=\frac{1}{n!}$. From this, the maximum permutation entropy $H_{\rm max}$ is quantified as

[TABLE]

Applying Eq. (2) for $n=3$ yields a maximum PE of approximately 2.585 bits. Using the maximum possible entropy $\log_{2}{n!}$, the normalized permutation entropy is calculated as

[TABLE]

Applying Eq. (3) to the example series $X$ results in $h_{3}\approx 0.5888$.

2 An overview of tools from TDA

Our TDA-based approaches for finding the delay dimension employs two types of persistence applied to two different types of data. Specifically, in the first approach we combine the [math]-D sublevel persistence of one-dimensional time series with the $z$-score, while in the second approach we utilize $1$-D persistent homology on embedded time series as part of the SW1PerS framework in Section 3.3. This section provides a basic background of the topics needed to at least intuitively understand the subsequent analysis. More specifics can be found in [17, 9, 13, 33].

2.1 Simplicial complexes

An abstract $k$-simplex $\sigma$ is defined as a set of $k+1$ indices where $\dim(\sigma)=k$. If we apply a geometric interpretation to a $k$-simplex, we can think of it as a set $V$ of $k+1$ vertices. Using this interpretation, a [math]-simplex is a point, a $1$-simplex is an edge, a $2$-simplex is a triangle, and higher dimensional versions can be similarly obtained.

A simplicial complex $K$ is a set of simplices $\sigma\subseteq V$ such that for every $\sigma\in K$, all the faces of $\sigma$, i.e., all the lower dimensional component simplices $\sigma^{\prime}\subset\sigma$ are also in $K$. For example, if a triangle ($2$-simplex) is in a simplicial complex $K$, then so are the edges of the triangle ($1$-simplices) as well as all the nodes in the triangle ([math]-simplices). The dimension of the resulting simplicial complex is given by the largest dimension of its simplices according to $\dim(K)=\max_{\sigma\in K}\dim(\sigma)$. The $n$-skeleton of a simplicial complex $K^{(n)}$ is the restriction of the latter to its simplices of degree at most $n$, i.e., $K^{(n)}=\{\sigma\in K\mid\dim(\sigma)\leq n\}$.

Given an undirected graph $G=(V,E)$ where $V$ are the vertices and $E$ are the edges, we can construct the clique (or flag) complex

[TABLE]

2.2 Homology

If we fix a simplicial complex $K$, then homology groups can be used to quantify the 1-dimensional topological features of the structure in different dimensions. For example, in dimension [math], the rank of the [math] dimensional homology group $H_{0}(K)$ is the number of connected components. The rank of the $1$-dimensional homology group $H_{1}(K)$ is the number of loops, while the rank of $H_{2}(K)$ is the number of voids, and so on. The homology groups are constructed using linear transformations termed boundary operators.

To describe the boundary operators, we let $\{\alpha_{\sigma}\}$ be coefficients in a field $F$ (in this paper we choose $F=\mathbb{Z}_{2}$). Since we are using the field $\mathbb{Z}_{2}$ we do not need to consider the orientation of the simplicial complex. Then $K^{(n)}$, the $n$-skeleton of $K$, can be used as a generating set of the $F$-vector space $\Delta_{n}(K)$. In this representation, any element of $\Delta_{n}(K)$ can be written as a linear combination $\sum\limits_{\sigma\in K^{(n)}}{\alpha_{\sigma}\sigma}$. Further, elements in $\Delta_{n}(K)$ are added by adding their coefficients. A finite formal sum of the $n$-simplices in $K$ is called an $n$-chain, and the group of all $n$-chains is the $n$th chain group $\Delta_{n}(K)$, which is a vector space.

Given a simplicial complex $K$, the boundary map $\partial_{n}:\Delta_{n}(K)\to\Delta_{n-1}(K)$ is defined by

[TABLE]

where $\hat{v}_{i}$ denotes the absence of element $v_{i}$ from the set. This linear transformation maps any $n$-simplex to the sum of its codimension $1$ (codim-$1$) faces. The geometric interpretation of the boundary operator is that it yields the orientation-preserved boundary of a chain.

By combining boundary operators, we obtain the chain complex

[TABLE]

where the composition of any two subsequent boundary operators is zero, i.e., $\partial_{n}\circ\partial_{n+1}=0$. An $n$-chain $\alpha\in\Delta_{n}(K)$ is a cycle if $\partial_{n}(\alpha)=0$; it is a boundary if there is an $n+1$-chain $\beta$ such that $\partial_{n+1}(\beta)=\alpha$. Define the kernel of the boundary map $\partial_{n}$ using $Z_{n}(K)=\{c\in\Delta_{n}(K):\partial_{n}c=0\}$, and the image of $\partial_{n+1}$ $B_{n}(K)=\{c\in\Delta_{n}(K):c=\partial_{n+1}c^{\prime},c^{\prime}\in\Delta_{n+1}(K)\}$. Consequently, we have $B_{k}(K)\subseteq Z_{k}(K)$. Therefore, we define the $n$th homology group of $K$ as the quotient group $H_{n}(K)=Z_{n}(K)/B_{n}(K)$. In this paper, we only need [math]- and $1$-dimensional persistent homology, and we always assume homology with $\mathbb{Z}_{2}$ coefficients which removes the need to keep track of orientation. In the case of [math]-dimensional homology, there is a unique class in $H_{0}(K)$ for each connected component of $K$. For $1$-dimensional homology, there is one homology class in $H_{1}(K)$ for each hole in the complex.

2.3 Filtration of a simplicial complex

Now, we are interested in studying the structure of a changing simplicial complex. We introduce a real-valued filtration function on the simplicies of $K$ such that $f(\tau)\leq f(\sigma)$ for all $\tau\leq\sigma$ simplices in $K$. If we let $\{y_{1}<y_{2}<\ldots<y_{\ell}\}$ be the set of the sorted range of $f$ for any $y\in\mathbb{R}$, then the filtration of $K$ with respect to $f$ is the ordered sequences of its subcomplexes

[TABLE]

The sublevel set of $K$ corresponding to $y$ is defined as

[TABLE]

where each of the resulting $K(y)$ is a simplicial complex, and for any $y_{1}\leq y_{2}$, we have $K(y_{1})\subseteq K(y_{2})$.

The filtration of $K$ enables the investigation of the topological space under multiple scales of the output value of the filtration function $f$. In this paper we consider two different filtration functions where each of these functions is applied to a different type of data: $1$-D persistence applied to point clouds embedding in $\mathbb{R}^{n}$ and [math]-D persistence applied to $1$-D time series.

$1$-D persistence applied to point clouds in $\mathbb{R}^{m}$:

A data type we use in this work is the set of vectors known as a point cloud $P\subset\mathbb{R}^{m}$. One common way to generate a point cloud $P$ is to embed a time series through the use of delay reconstruction. For each point $p\in P$, let $B(p,r)$ be the ball centered at $p$ and of radius $r$. Now, for some radius $r$, we can build a simplicial complex where, for example, the intersection of any two balls adds an edge, the intersection of three balls adds a triangle, and higher dimensional analogs are added similarly, see the construction for $r_{1}$ through $r_{4}$ in Fig. 3. The result of the union of all the balls $\cup_{p\in P}B(p,r)$ is called the Čech complex; although in practice we work with Rips complexes which are the clique complexes of the function $f$ and they are easier to obtain computationally while still capturing the interesting features of the space. The filtration function in this case is given by

[TABLE]

where $d(u,v)$ is the distance between vertices $u$ and $v$, and the distance between a vertex and itself is zero.

[math]-D persistence applied to $1$-D time series:

Let $\chi$ be the time-ordered set of the critical values of a time series. Here, we can think of the simplicial complex $K=G(V,E)$ containing a number of vertices $|V|$ equal to the number of critical values in the time series and only the edges $E$ that connect adjacent vertices. i.e., vertices $\{v_{i}\mid 1\leq i\leq n\}$, and edges $\{v_{i}v_{i+1}\mid 1\leq i\leq n-1\}$. Therefore, we have a one-to-one correspondence between the critical values in $\chi$ and the vertices of the simplicial complex $V$. We define the filtration function for every face $\sigma$ in $K$ according to

[TABLE]

Using this filtration function in Eq. (4), we can define an ordered sequence of subcomplexes where $y\in[\min{(\chi)},\max{(\chi)}]$.

2.4 Persistent homology

Persistent homology is a tool from topological data analysis which can be used to quantify the shape of data. The main idea behind persistent homology is to watch how the homology changes over the course of a given filtration.

Fix a dimension $n$, then for a given filtration

[TABLE]

we have a sequence of maps on the homology

[TABLE]

We say that a class $[\alpha]\in H_{n}(K_{i})$ is born at $i$ if it is not in the image of the map $H_{n}(K_{i-1})\to H_{n}(K_{i})$. The same class dies at $j$ if $[\alpha]\neq 0$ in $H_{n}(K_{j-1})$ but $[\alpha]=0$ in $H_{n}(K_{j})$.

This information can be used to construct a persistence diagram as follows. A class that is born at $i$ and dies at $j$ is represented by a point in $\mathbb{R}^{2}$ at $(i,j)$. The collection of the points in the persistence diagram, therefore, give a summary of the topological features that persists over the defined filtration. See the example of Fig. 4 for $n=0$ and time series data, and Fig. 3 for $n=1$ and point cloud data.

3 Permutation Delay

To form permutations from a time series a delay embedding is applied to a uniform subsampling of the original time series according to the embedding parameter $\tau$. For example, the subsampled sequence $X$ with elements $\{x_{i}:i\in\mathbb{N}\cup 0\}$ subject to the delay $\tau$ is defined as $X(\tau)=[x_{0},x_{\tau},x_{2\tau},\ldots]$. Riedl et al. [40] showed that PE is sensitive to the time delay, which prompts the need for a robust method for determining an appropriate value for $\tau$. For estimating the optimal $\tau$, we will be investigating the following methods in the subsequent sections: Mutual Information (MI) in Section 3.1, combining [math]-D persistence with the frequency and time domains (Section 3.4), and SW1PerS (Section 3.3). We recognize, but do not investigate, some other commonly used methods for finding $\tau$. These include the autocorrelation function [18] and the phase space expansion [8].

Before introducing the TDA-based methods, in Section 3.2 we elucidate a connection between the permutation delay parameter and the state space reconstruction delay parameter for Takens’ embedding to justify the use of methods such as MI and the TDA-based methods.

3.1 Mutual Information

A common method for selecting the delay $\tau$ for state space reconstruction is through a measure of the mutual information between a time series $x(t)$ and its delayed version $x(t+\tau)$ [15]. Mutual information is a measurement of how much information is shared between two sequences, and it was first realized by Shannon et al. [43] as

[TABLE]

where $X$ and $Y$ are separate sequences, $p(x)$ and $p(y)$ are the probability of the element $x$ and $y$ separately, and $p(x,y)$ is the joint probability of $x$ and $y$. Fraser and Swinney [15] showed that for a chaotic time series the mutual information between the sequence $x(t)$ and $x(t+\tau)$ will decrease as $\tau$ increases until reaching a minimum. At this delay $\tau$, the individual data points share minimum amount of information, thus indicating that the data points are sufficiently separated. While this delay value was specifically developed for phase space reconstruction from a single time series, we show in Section 3.2 how this delay parameter selection method is also appropriate for permutation entropy. Due to this relationship and the result in [31] suggesting the use of MI for PE delay selection, we benchmark the delays $\tau$ obtained from our methods against their MI counterparts.

3.2 Relating Permutations to Delay Reconstruction

The goal in this section is to relate the distribution of permutations formed from a given delay $\tau$ to the state space reconstruction with the same delay $\tau$. This connections will show the time delay for both permutations and state space reconstruction are related.

Let us first start by describing the process for state space reconstruction and its similarity to permutations. As described by Takens’ [46], we can reconstruct an attractor that is topologically equivalent to the original attractor of a dynamical system through embedding a 1-D signal into $\mathbb{R}^{n}$ by forming a cloud of delayed vectors as $v_{i}=[x(t_{i}),x(t_{i+\tau}),x(t_{i+2\tau}),\dots,x(t_{i+(n-1)\tau})]$ for $i\in[0,L-n\tau]$, where $L$ is the length of the discretely and uniformly sampled signal. Permutations are formed in a very similar fashion where we take our vectors $v_{i}$ and find their symbolic representation based on their ordinal ranking as explained in Section 1.1. The different permutation types can be viewed as a inequality-based binning of the $\mathbb{R}^{n}$ vector space of the reconstructed dynamics as shown in Fig. 5 for dimension $n=3$. This provides a first intuitive understanding of the connection between permutation and state space reconstruction; however, we need to determine if the optimal $\tau$ parameter used in $v_{i}$ is related to the optimal delay in PE.

Takens’ embedding theorem explains that, theoretically, any delay $\tau$ would be suitable for reconstructing the original topology of the attractor; however, this has the requirement of unrestricted signal length with no additive noise in the signal [46]. Since this is rarely a condition found in real-world signals, a $\tau$ is chosen to unfold the attractor such that noise has a minimal effect on the topology of the reconstructed dynamics.

Let us now explain what we mean by the correspondence between $\tau$ and the unfolding of the dynamics and what effect this has on the corresponding permutations. If the delay $\tau$ is too small (e.g. $\tau=1$ for a continuous dynamical system with a high sampling rate) the delay embedded reconstructed attractor will be clustered around the hyper-diagonal in $\mathbb{R}^{n}$. Additionally, the corresponding permutations will be overwhelmingly dominated by the permutation types $\pi_{1}$ and $\pi_{n!}$ with these two permutations being of the all increasing and all decreasing ordinal patterns, respectively. The dominance of these two permutations for a delay $\tau$ that is too small was termed by De Micco et al [11] as the “redundancy effect." For an example of this see the permutation distribution and clustering about the hyper-diagonal in $\mathbb{R}^{3}$ in Fig. 6 when $\tau=1$. This example is based on the $x$-solution to the periodic Rossler dynamical system as described in Section A.2. As the delay increases beyond the redundancy effect, the reconstructed attractor begins to unfold to have a similar shape and topology as the true attractor. Correspondingly, as the delay increases the permutation distribution tends towards a more equiprobable distribution (See Fig. 6 at $\tau\approx 14$). A way of summarizing the permutation probability distribution is actually through PE itself and more specifically the analysis of Multi-scale Permutation Entropy (MsPE). Riedl et al. [40] showed how after the redundancy effect there is a suitable delay for PE, which we related to the first maxima of the MsPE plot [31]. The MsPE plot for our periodic Rossler example is shown in Fig. 6. Let us also look at the MI plot as a comparison. The idea behind MI is that at the first minima of the mutual information between $x(t_{i})$ and $x(t_{i+\tau})$ the delay $\tau$ accurately provides a suitable delay for state space reconstruction. A quick investigation of the MI function reveals a high degree of correlation between the MI function and the MsPE function with the first maxima of MsPE being approximately at the same $\tau$ as the first minima of MI. When the delay becomes significantly larger than the first minima of MI or maxima of MsPE, the permutation distribution begins to fluctuate as shown in Fig. 6. This effect was termed the “irrelevance effect" by De Micco et al. [11]. This increasing of $\tau$ beyond the the first minima also correlates with what Kantz and Schreiber [22] describe as the reconstruction filling an overly large space with the vectors already being independent. Additionally, at a minima beyond the first minima, Fraser and Swinney [15] showed how the reconstructed attractor shape will no longer qualitatively match the shape of the true state space.

In summary, we have shown a main heuristic result we need to move forward: tools for delay parameter selection for PE can be suitable for state space reconstruction and vice-versa. While we do not provide a proof that PE and state space reconstruction use the same $\tau$, it has recently been shown that there is a connection between co-homology, information theory, and probability does exist [7], which strengthens our qualitative analysis of this connection. In the following sections we leverage tools from TDA to determine the optimal $\tau$ associated with the unfolding of the attractor.

3.3 Finding $\tau$ Using Persistent Homology

In this section we develop a novel method based on persistent homology for estimating an appropriate delay for permutations and state space reconstruction. Specifically, we investigate the effects of varying $\tau\in[1,\tau_{\max}]$ on the calculation of the maximum persistence and the periodicity score from SW1PerS [36]. Perea and Harer developed SW1PerS as a TDA method for measuring periodicity in a time series; however, our goal is to leverage this method to determine a suitable $\tau$ for permutation entropy and state space reconstruction based on the unfolding of an attractor and the associated 1-D persistent homology.

SW1PerS uses 1-D persistent homology as described Section 2 to measure how periodic or significant the circular shape of an embedded time series (point cloud) is as $\tau$ increases, which corresponds to the embedding window size increasing as $w=m\tau$ with $m\in\mathbb{N}$ as the embedding dimension of the sliding window vector. Specifically, the sliding window $SW$ for SW1perS is defined as

[TABLE]

where $f(N)(t)$ is a truncated Fourier series of the signal to $N$ Fourier series terms and $\tau$ and $m$ are, respectively, SW1PerS’ embedding delay and dimension. Applying Eq. (6) to a sliding window of width $w$ across the domain of the time series results in a collection of vectors known as a point cloud, which live in an $m$-dimensional Euclidean space. However, it may not be desirable to use all of the embedded vectors from (6) due to the $O(n^{3})$ time complexity of calculating the persistent homology of a point cloud via the Vietoris-Rips complex. To improve the calculation time we chose to use a sparse version of the point cloud through a subsampling to have $n_{T}$ windows from the original point cloud. We chose to set the number of sliding windows as $n_{T}=200$ to be sufficiently high to detect circular structure in the embedding.

For SW1PerS, Perea et al. [36] determined $m$ using Shannon-Nyquist sampling theoem according to $m\geq 2N$ (here we use $m=2N$), where N is the number of Fourier series terms which was arbitrarily chosen in [36]. In contrast, we automate choosing $N$ by approximating the Fourier series using the discrete Fourier transform. To do this we compute the normalized $\ell_{2}$ norm between the reconstructed signal from the truncated Fourier series and the original signal. The $\ell_{2}$ norm is used to obtain the value of $N$ that yields an error within a desired threshold of $\ell_{2}(N)<0.25$. Specifically, if we let the time series $X$ be a discrete time sampling of a piece-wise smooth signal $x(t)$, then the $N$-partial sum of the Fourier series of $x(t)$ can be approximated according to

[TABLE]

where $X$ is the original signal that has been point-wise centered and normalized with $|X|$ as the length of the signal. As a rule of thumb $N\approx|X|/8$ yields an accurate reconstruction of $x(t)$ [4], which we use as an upper bound of $N$. The relative $\ell_{2}$ norm that measures the error between time series $X$ and its reconstruction $f(N)(t)$ is given by

[TABLE]

For our application, we consider $f(N)(t)$ as sufficiently close to $x(t)$ when we find a value of $N$ for which $\ell_{2}(N)<0.25$. We chose $0.25$ because it leads to embedding dimensions $m$ which are not overly large ($m<10$ typically), and it deals with the possibility of moderate additive noise in the signal. Using the truncated Fourier series we are also able to determine an upper bound for $\tau$ using the Nyquist sampling criteria as

[TABLE]

where $f_{\rm sig}$ are the $N$ significant frequencies from the truncated fast Fourier transform. This gives in an automatic way all the components we need to apply SW1PerS.

To determine the optimal delay using persistent homology, we investigated two summaries of the resulting persistence diagrams from SW1PerS: (1) the maximum lifetime

[TABLE]

with $\tilde{\mathcal{D}}$ as the SW1PerS persistence diagram and (2) the periodicity score, which was defined in [35] as

[TABLE]

where $r_{B}$ and $r_{D}$ are the birth and deaths times associated to $\max({\rm pers}(\tilde{\mathcal{D}}_{1})$. We then calculate these point summaries over the range $\tau\in[1,\tau_{\max}]$ to generate $\vec{s}$ and $\vec{\mathcal{L}}$ for the periodicity scores $s$ and persistence maxima $\mathcal{L}$, respectively.

To demonstrate the functionality of this method, consider a simple example using the periodic Rossler system as described in Eq. (23) (see Fig. 7). This example shows three different window sizes for embedding dimension $m=2$ (this dimension was chosen for visualization purposes) and $\tau=1$, 10, and 17, to show the resulting scores for small, optimal and overly large window sizes, respectively. Figure 7 shows that the optimal window size at $\tau_{o}=10$ results in a maximum $\mathcal{L}$ and minimum $s$ over the range $\tau\in[1,\tau_{\max}]$, where $\tau_{\max}=20$ from the truncated Fourier spectrum. This suggests that an appropriate delay for both state space reconstruction via Takens’ embedding and permutation entropy is $\tau=10$.

For a chaotic time series, choosing $\tau$ from the minimum or maximum of $s$ and $\mathcal{L}$ is not as trivial as the example shown in Fig. 7. Specifically, due to the non-linear behavior of a chaotic time series there may not always be a clear, single minima as shown in the example periodic Rossler system, but rather two or more local minima with similar prominence. To accurately approximate the average minima and select an associated delay $\tau$, we will use heuristic cutoffs $C_{s}$ and $C_{\mathcal{L}}$, where these cutoffs are defined as $C_{s}={\frac{1}{2}[\max(s)+\min(s)]}$ and $C_{\mathcal{L}}={\frac{1}{2}[\max(\mathcal{L})+\min(\mathcal{L})]}$. Specifically, we will choose $\tau$ based on the average $\tau$ such that $s\geq C_{s}$ or $\mathcal{L}\geq C_{\mathcal{L}}$. To demonstrate this method we use a chaotic response of the Rossler system from Eq. (23) and calculate the two cutoffs as shown in Fig. 8. This example results in an average delay greater than $C_{\mathcal{L}}$ as $\tau=12$ and less than $C_{s}$ as $\tau=12$. This example demonstrated that the method of selecting the average $\tau$ greater or less than the cutoffs results in a similar $\tau$ for both periodic and chaotic time series.

3.4 Finding $\tau$ Using Sublevel Set Persistence

In this section our goal will be to leverage sublevel set persistence for the selection of $\tau$ for both state space reconstruction and permutation entropy. Specifically, our goal is to automate the frequency analysis method [29] for selecting $\tau$ for state space reconstruction by analyzing both the time and frequency domain of the signal using sublevel set persistence. Melosik and Marszalek [29] leveraged Shannon-Nyquist sampling criteria and used the maximum significant frequency $f_{\rm max}$ and the sampling frequency $f_{s}$ to select an appropriate $\tau$ according to

[TABLE]

where $\alpha\in[2,4]$. The value $\alpha=2$ is associated to the Nyquist sampling rate, while $\alpha>4$ produces an oversampling. Since this method was developed using the Nyquist sampling rate, it is applicable continuous, band-limited signals. This frequency based approach was used to find suitable delays for the 0/1 test on chaos and heuristically compare the Lorenz attractor and its time-delay reconstruction [29]. The heuristic comparison showed that this frequency approach actually provided more accurate delay parameter selections for state space reconstruction than the mutual information function when trying to replicate the shape of the attractor. Unfortunately, a major drawback of this method is the non-trivial selection of $f_{\max}$. In Melosik’s and Marszalek’s original work [29] the maximum frequency was manually selected using a normalized, such that $\in[0,1]$, Fast Fourier Spectrum (FFT) cutoff of approximately 0.01, which does not address the possibility of additive noise.

In our previous work [31] we approximated the maximum “significant" frequency in a time series using the FFT and defining a power spectrum cutoff based on the statistics of additive noise in the FFT. An issue with this method for non-linear time series is that the Fourier spectrum does not easily yield itself to selecting the maximum “significant" frequency for chaotic time series even with an appropriately selected cutoff to ignore additive noise. Additionally, the method was only developed for Gaussian White Noise (GWN) contamination of the original time series.

The following sections improve the selection of the maximmum significant frequency using two novel methods based on 0-D sublevel set persistence. We chose to use 0-D sublevel set persistence due to its computational efficiency and stability for true peak selection [24, 10]. The first method is based on a time domain analysis of the sublevel set lifetimes (see Section 3.4.1), and the second implements a frequency domain analysis using sublevel set persistence and the modified $z$-score (see Section 3.4.2).

3.4.1 Time Domain Approach

The first approach we implement for estimating the maximum significant frequency of a signal is based on a time domain analysis of the sublevel set persistence. This process uses the time ordered lifetimes from the sublevel set persistence diagram. We previously introduced time ordered lifetimes and a cutoff separating the sublevel sets associated with noise in [5]. Here we use those methods and results to find the time $t_{B}$ in which all the significant sublevel sets are born. Figure 9 shows an example time order lifetimes plot where the time between two adjacent lifetimes is defined as $T_{B_{i}}$.

If we use $T_{B_{i}}$ as an approximation of a period in the time series, then we can calculate the associated frequencies as $f_{i}=1/T_{B_{i}}$ Hz. If we then look at the distribution of $f_{i}$, the maximum “significant" frequency can be approximated using the 75% quantile of the distribution of the frequencies as $f_{\rm max}\approx Q_{75}(f)$. This quantile allows for a few outlier frequencies to occur without having a significant effect on the estimate of the maximum significant frequency.

Applying this method to the periodic Rossler system described in Eq. (23) results in $\tau=10$ with the corresponding state space reconstruction for $n=2$ shown in Fig. 10. This suggested delay is close to that of mutual information ($\tau=12$). This suggests that the time-domain analysis for selecting the maximum frequency and corresponding delay functions can automatically and accurately suggest an appropriate delay for permutation entropy and state space reconstruction.

3.4.2 Fourier Spectrum Approach

In this section we present a novel TDA based approach for finding the noise floor in the Fourier spectrum for selecting the maximum significant frequency $f_{\max}$ to be used for selecting $\tau$ for PE through Eq. (12). Specifically, we show how the [math]-dimensional sublevel set persistence, a tool from TDA discussed in Section 2, can be used to find the significant lifetimes and the associated frequencies in the frequency spectrum. Although it would be ideal to separate the significant lifetimes based on propagating the FFT of a random process into the persistence space, this task is not trivial. There have been studies on pushing forward probability distributions into the persistence domain [1, 2, 21], but it is difficult to obtain a theoretical cutoff value in persistence space; therefore, we instead seperate the noise lifetimes from significant lifetimes through the use of the modified $z$-score. This separation allows us to find the noise floor and the maximum significant frequency via a cutoff. This process for finding the cutoff and associated maximum frequency is illustrated in Fig. 11. The following paragraphs give an overview of the modified $z$-score and cutoff analysis.

Modified $z$-score

The modified $z$-score $z_{m}$ is essential to understanding the techniques used for isolating noise from a signal [41]. The standard score, commonly known as the $z$-score, uses the mean and the standard deviation of a dataset to find an associated $z$-score for each data point and is defined as

[TABLE]

where $x$ is a data point, $\mu$ is the mean, and $\sigma$ is the standard deviation of the dataset, respectively. The $z$-score value is commonly used to identify outliers in the dataset by rejecting points that are above a set threshold, which is set in terms of how many standard deviations away from the mean are acceptable. Unfortunately, the $z$-score is susceptible to outliers itself with both the mean and the standard deviation not being robust against outliers [25]. This led Hampel [19] to develop the modified $z$-score as an outlier detection method that is robust to outliers. The logic behind the modified $z$-score or median absolute deviation (MAD) method is grounded on the use of the median instead of the mean. The MAD is calculated as

[TABLE]

where $\mathbf{x}$ is an array of data points and $\tilde{x}$ is its median. The MAD is substituted for the standard deviation in Eq. (13). To improve the modified $z$-score, Iglewicz and Hoaglin [20] suggested to additionally substitute the mean with the median. The resulting equation for the modified $z$-score is then

[TABLE]

where the value 0.6745 was suggested in [20]. We can now use the modified $z$-score $z_{m}$ for evaluating the “significance" of each point in the sublevel set persistence diagram of the Fourier spectrum. A threshold for separating noise in the persistence domain is discussed in the following paragraph.

Threshold and Cutoff Analysis

To determine the noise floor in the normalized Fast Fourier Transform (FFT) spectrum, we compute the [math]-dimensional persistence of the FFT. This provides relatively short lifetimes for the noise, while the prominent peaks, which represent the actual signal, have comparatively long lifetimes or high persistence. To separate the noise from the outliers we calculate the modified $z$-score for the lifetimes in the persistence diagram. We can then determine if the lifetime is associated to noise or signal based on a $z_{m}$ cutoff as $D$, where we can label a lifetime as signficant (an outlier) if $z_{m}>D$. Iglewicz and Hoaglin [20] suggest a $z_{m}$ threshold of $D=3.5$ based on an analysis of 10,000 random-normal observations. However, we apply both the FFT and 0-D sublevel set persistence to the original signal so we need to determine if this cutoff is suitable for our application. To do this we used a signal of 10,000 random-normal observations and applied FFT. We then calculated the 0-D sublevel set persistence and computed the modified $s$-score $z_{m}$ using the resulting lifetimes. For an accurate cutoff we would expect to label all of the lifetimes as noise with $z_{m}<D$ since each signal is composed of pure noise. As shown in Fig. 12, a threshold of approximately $D=4.8$ labels all of the lifetimes as noise. This threshold was rounded up to 5 for simplicity. We can now simply define a cutoff based on the labeling of each lifetime from the modified $z$-score with $\rm Cutoff=\max(lifetime_{noise})$.

We can now find the maximum significant frequency $f_{\rm max}$ as the highest frequency in the Fourier spectrum with an amplitude greater than the specified cutoff. For this method to accurately function, it is required that there is some additive noise in the time series. To accommodate this, additive Gaussian noise with Signal-to-Noise Ratio of 30 dB is added to the time series before calculating the FFT. If we apply this method to the example periodic Rossler system time series we find a suggested delay of $\tau=5$. In comparison to mutual information this delay is approximately half as large as it should be. However, we will further investigate its accuracy on several other systems in Section 5 to make more general conclusions on the functionality of this method for selecting $\tau$.

4 Permutation Dimension

In this section we will show that, contrary to the delay selection, the dimension for permutation entropy is not related to that of Takens’ embedding. Additionally, we will provide a simple method for selecting an appropriate permutation dimension based on the permutation distribution.

Permutation entropy is often used to differentiate between the complexity of a time series when there is a dynamic state change (e.g. periodic compared to chaotic), so the dimension should be chosen such that it is large enough to capture these changes. To accomplish this we suggest that permutations of the time series do not occupy all of the possible permutations, but rather only a fraction of the permutations when an appropriate delay is selected. This criteria is set so that a change can be captured by an increase/decrease in the number of permutations and their associated probabilities. Because of this, we suggest a dimension where, at most, only 50% of the permutations are used. However, it may be more suitable to select a dimension where a lower percent is used (e.g. 10%).

To begin this method for determining if the dimension is high enough to capture the time series complexity we will define $N_{\pi}$ as the number of permutation types where the probability of that permutation type is significant. Specifically, we will consider the probability of that permutation to be significant if the number of occurrences of permutation $\pi$ is greater than 10% of the maximum number of occurrences of any permutation type of dimension $n$. The permutation delay $\tau$ was selected from the expert suggested values provided in [31, 40]. We can now express our needed dimension as the ratio and inequality

[TABLE]

where $R=0.50$ for the suggested maximum 50% criteria.

To compare this dimension to the standard Takens’ embedding tools for selecting $n$ we will implement four examples:

[TABLE]

where $t\in[0,100]$ with a sampling rate of 20 Hz and $\mathcal{N}$ is Gaussian additive noise. By applying Eq. (16) to the time series in Eq. (17), we can suggest dimensions of 2 for $x_{1}(t)$, 4 for $x_{2}(t)$, 6 for $x_{3}(t)$, and 7 for $x_{4}(t)$ as shown in Fig. 13.

In comparison to Takens’ embedding, for time series $x_{2}(t)$ dimension $n=2$ would be sufficient, but if this was used for permutation entropy, no increase in complexity could be detected. Additionally, this result suggests an upper bound on the dimension for permutation entropy as $n\approx 9$ as the ratio in Eq. (16) is approximately 0 for dimensions $n>9$. As a rule of thumb from this result, a dimension of $8$ would be suitable for almost all applications, but it would be optimal to minimize the dimension to reduce the computation time of PE. In Section 5 we will show the resulting suggested dimensions using this method for a wide variety of dynamical systems.

5 Results

This section provides the results of the parameter selection methods. First, in Section 5.1, we calculate the delay parameter for a wide variety of dynamical systems and data sets using mutual information and the the automatic TDA-based methods described in this manuscript. Unfortunately, the optimal parameters cannot be decided based on a simple entropy value comparison since there is no direct equivalence between PE and other entropy approximations of a signal such as Kolmogorov-Sinai (KS) entropy with only a bounding between the two as $\rm KS\leq PE$ [23]. Therefore, to determine the accuracy of the automatically selected PE parameters we implement two other methods of comparison. The first is a comparison to expert suggested parameters for a wide variety of systems (see Section 5.1). The second approach is a comparison to optimal parameters based on having a significant difference between the PE of two different states for each system. Of course the second method has the requirement that we have a system model or data set with two different states for comparison, which is not typically the case, but does allow for an approximation of optimal PE parameters for these systems. These comparisons are discussed in Section 5.1.

The second half of the results, in Sections 5.2 and 5.3, is based on analyzing the robustness of the automatic TDA-based PE parameter selection methods to additive noise and signal length, respectively.

5.1 Parameter Value Comparison for Common Dynamical systems

To determine a range of approximately optimal PE parameters we will quantify the difference between PE values for a wide range of delays and dimensions with the difference for a given $\tau$ and $n$ calculated as

[TABLE]

where the superscripts Ch. and Pe. represent the PE calculation on the chaotic and periodic time series for the given dynamical system. The specific parameters used to generate periodic and chaotic responses for each system are described in Appendix A. If we apply Eq. (18) to the Rossler system for $\tau\in[1,15]$ and $n\in[3,10]$ we find that $\Delta h_{n}(\tau)$ is significant when $\tau\in[9,15]$ and $n\in[6,10]$ as shown in Fig. 14. However, as mentioned previously in Section 4, dimensions greater than 8 can be computationally expensive.

We consider this range where $\Delta h_{n}(\tau)$ is relatively large as the range of optimal PE parameters to be compared to. We repeated this process for finding the optimal parameter ranges for PE using a similar procedure to this Rossler example as shown in Table 1.

To verify our TDA-based methods for determining $\tau$, Table 1 compares our results to the values from a wide variety of systems for both the first minima of the mutual information function and from expert suggestions, including several listed by Riedl et al. [40]. The table also shows the resulting permutation dimensions suggested from the permutation statistics as described in Section 4 for both $R=0.1$ and $R=0.5$ from Eq. (16). For these systems we have also included, where applicable, the delay and dimension parameter estimates for both periodic and chaotic responses to validate each method’s robustness to chaos and non-linearity. However, for the medical data section we instead included a healthy/control and unhealthy (arrhythmia for ECG and seizure for EEG) as a substitute for a periodic and chaotic response, respectively. A detailed description of each dynamical system or data set used, including parameters for periodic and chaotic responses, is provided in the Appendix.

In table 1 we have highlighted the methods that failed to provide an accurate delay $\tau$ in red. We will now go through the methods and highlight the advantages and drawbacks as well as provide general suggestions for which method to use based on the category.

Noise Models: We only have one expert suggestion of parameters for the noise models category, which is for Gaussian white noise (Gauss.) as $\tau=1$ and $n\in[3,6]$. In regards to the delay, all TDA based methods show an accurate selection of $\tau=1$, however the suggestion of $\tau=3$ from Mutual Information (MI) is slightly higher than suggested. We found that the expert suggested dimensions of 3 to 6 is significantly lower than the minimum dimension suggested by our permutation statistics method of $n=7$. As mentioned in Section 4, we believe it is necessary to have the number of permutation used to be atleast less than 50% of all the permutation available, which corresponds to a dimension $n=7$ for Gaussian noise. Additionally, as shown in Fig. 13, if a dimension of even $n=6$ is used approximately 95% of the permutations are already in use for Gaussian noise, which would make dynamic state changes difficult or statistically insignificant if there is a complexity increase. From this logic we can then conclude that a suitable dimension should actually be atleast $n=7$ if any increase in the time series complexity is expected. If only decreases in complexity are expected, then a dimension of $n=6$ may be suitable.

Continuous Flows: The next category is of continuous flows described by systems of non-linear differential equations. As shown in Table 1, both the time domain analysis via sublevel set persistence and mutual information provide accurate delay suggestions for all of the examples. However, the 1-D persistent homology methods discussed in Section 3.3 also provide an accurate delay for every systems besides for the chaotic Chua circuit. This shortcoming was most likely due to an inaccurate selection of the maximum significant frequency and associated $\tau_{\max}$. We can also conclude that the frequency domain analysis using sublevel set persistence consistently provided delays that were too small. In regards to the dimension, the suggested dimensions from the permutation statistics agreed with the delay suggested by experts for all of the continuous flow systems. This suggests that our method for selecting a dimension for permutation entropy in Section 4 is accurate for simulations of continuous differential equations.

Periodic Functions: For periodic functions, including a simple sinuisodal function (periodic) and two incommensurate sinuisoidal functions (quasiperiodic), our results in Table 1 show that all methods, including mutual information, provide accurate selections of $\tau$ except the Fourier spectrum analysis via sublevel sets. This method results in a significantly high suggestion for $\tau$. In regards to the dimension selection, our results using the permutation statistics method described in Section 4 agree with the expert suggested minimum dimension of $n=4$.

Maps: When selecting the delay parameter for permutations and takens’ embedding for maps we found that all of the topological methods suggested accurate delay parameters, while the standard mutual information method selected overly large delay parameters when the maps are chaotic. Therefore, we suggest the use of one of the topological methods when estimating the delay parameter for maps. For the permutation dimension we found a suggested dimension $n\in[4,7]$, in comparison to the expected suggested dimension ranging from 2 to 16. While the range suggested from the permutations statistics as described in Section 4 falls within the range suggested by experts, their range is too broad. Specifically, a dimension greater than 9 can be computationally cumbersome, and a dimension lower than 4 would not show significant differences for dynamic state changes. Therefore, we suggest the use of our narrower range of dimension from $n\in[5,6]$ for maps which agrees with our optimal PE parameter range.

Medical Data: The medical data used in this study Inherently has some degree of additive noise, which provides a first glimpse into the noise robustness of the delay parameter selection methods investigated. However, a more thorough investigation will be provided in Section  5.2. From our analysis, we disagree with the expert suggested delay $\tau\in[1,3]$, but rather the delay selected from either mutual information or the time domain analysis of sublevel set persistence. The general selection for delays between 1 and 3 does not account for the large variation in possible sampling rates. If a small delay is used in conjunction with a high sampling rate, an inaccurate delay could be selected resulting in indistinguishable permutation entropy values as the dynamic state changes. In regards, to the permutation dimension $n$, we believe that a more appropriate dimension, in comparison to the values suggested by experts, should range between 5 and 7 for medical data applications.

5.2 Robustness to Additive Noise

To determine the noise robustness of the delay parameter selection methods investigated in this work we will use an example time series. Specifically, we will use the $x$ solution to the periodic Rossler system as described in Section A.2. We will use additive Gasussian noise $\mathcal{N}(\mu=0,\sigma^{2})$, where $\sigma$ is determined from the Signal-to-Noise Ratio (SNR). The SNR is a measurement of how much noise there is in the signal with units of decibels (dB) and is calculated as

[TABLE]

where $A_{\rm signal}$ and $A_{\rm noise}$ are the Root-Mean-Square (RMS) amplitudes of the signal and additive noise, respectively. If we manipulate Eq. (19) we can solve for $A_{\rm noise}$ as

[TABLE]

Because $x(t)$ is a discrete sampling from a continous system with $t=[t_{1},t_{2},\ldots,t_{N}]$, we calculate $A_{\rm signal}$ as

[TABLE]

where $\bar{x}$ is the mean of $x$ and is subtracted from $x(t)$ to center the signal about zero. with $A_{\rm noise}$ calculated, we set the additive noise standard deviation as $\sigma=A_{\rm noise}$.

We applied a sweep of the SNR from 1 to 40 in increments of 1 with each SNR being repeated for 30 unique realizations of the noise distributed as $\epsilon\sim\mathcal{N}(0,A_{\rm noise}^{2})$. For each realization of $x(t)+\epsilon$ the delay parameters were calculated using all 5 methods: sublevel set persistence of the frequency domain $\tau_{\rm SL_{f}}$, sublevel set persistence of the time domain $\tau_{\rm SL_{t}}$, the minima of SW1PerS score $\tau_{\rm PH_{s}}$, the maxima of the maximum persistence $\tau_{\rm PH_{\mathcal{L}}}$, and mutual information $\tau_{\rm MI}$. The mean and standard deviation of the 30 trials at each SNR were calulcated for each method as shown in Fig. 15.

Figure 15 shows that the sublevel set persistence methods fail to provide an accurate delay $\tau$ in comparison to the expert suggested delay $\tau_{\rm exp.}=9$ when $\rm SNR<10$ dB. While this does show a limit for the sublevel set persistence methods, SNR values below 10 dB indicte extremely noisy signals. In signal processing, a rule-of-thumb is to consider the signal to contain extractable information if SNR > 15 dB. However, the 1-D Persistent homology methods and mutual information provide accurate delay parameter selection down to an SNR of 2 dB.

5.3 Robustness to Signal Length

A common issue with signal processing and time series analysis methods is their limited functionality with smaller sets of data available, which has been used to analyze the sentitivty of the delay parameter selection [12]. Here we will investigate the limitations of these methods in the face of short time series. We will do this analysis by incrementing the length of the time series with the PE parameters calculated at each increment. For our analysis we will again use the Rossler system as described in Section A.2. Specifically, we increment the length of the signal from $L=75$ to 1000 in steps of 25 (see Fig. 16. However, if this type of analysis is not available for the data set being analyzed, for time series analysis applications it is commonly suggested to have a data length of $L=4000$ for continuous dynamical systems and and $L=500$ for maps [49].

In Fig. 16 we see that all of the methods reach an accurate value of $\tau$, in comparison to the expert suggested $\tau=9$ when the time series contains at least 125 data points. An important note to make is that this result is not general for all continuous dynamical systems. The required length of the signal is going to vary depending on the sampling rate of the time series. To determine a general requirement for the methods we repeated this analysis method for all of the systems shown in Table 1. Our result from this analysis found that, in general, a signal of $L\geq 15\tau$ allows selecting an appropriate PE and state space reconstruction delay $\tau$ using the TDA-based methods described in this manuscript.

6 Conclusion

We described two novel TDA-based approaches for automatically determining the PE delay parameter $\tau$ given a sufficiently sampled/oversampled time series. The first TDA method is based on SW1PerS, which operates on a sliding window embedding of the time-domain signal. By varying $\tau$ and calculating the persistent homology of the sliding window embedding at each $\tau$ we were able to show how two different statistical summaries of the persistence diagram could be used to select an appropriate delay for PE. Specifically, we implemented the periodicity score from SW1PerS, and the maximum persistence, where an appropriate delay is associated with either the first maximum of the maximum persistence or the first significant minima of the periodicity score. The other TDA based methods implemented sublevel set persistence. We investigated the sublevel set persistence in both the time and frequency domains to determine the maximum signficant freuqnecy, which was then used to estimate an appropriate delay based on the Shannon-Nyquist sampling criteria [29].

In regards to the permutation dimension $n$, in Section 4 we developed a simple statistical analysis method for selecting an appropriate permutation dimension $n$ based on the need for only a portion of the permutations to be used in the time series to capture complexity changes. This method also revealed that the permutation dimension and Takens’ embedding dimension are not necessarily related and tools for the Takens’ embedding dimension cannot generally be used for the permutation dimension.

To determine the accuracy of these methods, the resulting delays were compared to each of the standard MI approach of Fraser and Swinney [15], expert suggested parameters of various categories of dynamical systems and data sets, and optimal parameters calculated if the dynamical system model allowed for various differing dynamical states to be simulated. This result showed that the sublevel set persistence of the time domain method provided the most accurate delays for all of the systems analyzed except for EEG data. However, due to the differing sampling rates for different EEG data sets, the expert suggested delay of 1–3 may not be accurate for the EEG data investigated in this work. As an alternative method, the SW1PerS approach for delay parameter selection provided accurate delay suggestions except for the chaotic Chua circuit, which had delays that were significantly too large. This may be due to the need for a higher dimension analysis of the sliding window embedding as only the 1-D features (loops) were investigated. We do not recommend using the sublevel set persistence of the frequency domain as this method consistently provided delays that were too small for continuous flow differential equations and too large for periodic functions. Our noise robustness analysis revealed that all of the methods we developed here were robust to additive noise up to the SNR values as low as 10 dB. We also analyzed the robustness of the methods to signal length and found they provide an accurate delay even for short time series with the general suggestion of signal length $L>15\tau$.

In comparison to the expert suggested dimensions and the optimal dimensions from comparing chaotic and periodic signals, our dimension parameter selection method accurately provided similar for all of the systems investigated in this work. Additionally, the range suggested using our method was more precise than the dimensions suggested by experts giving the user a more definite answer to an appropriate permutation dimension.

Acknowledgment

This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-22-1-0007.

Appendix A Summary of the used data and models

A.1 Lorenz System

The Lorenz system used is defined as

[TABLE]

The Lorenz system had a sampling rate of 100 Hz. This system was solved for 100 seconds and the last 20 seconds were used. For a periodic response parameters $\sigma=10.0$, $\beta=8.0/3.0$, and $\rho=100$ were used. For a chaotic response parameters $\sigma=10.0$, $\beta=8.0/3.0$, and $\rho=105$ were used.

A.2 Rössler System

The Rössler system used was defined as

[TABLE]

with parameters of $a=0.10$ for periodic and $a=0.15$ for chaotic, $b=0.20$, $c=14$, which was solved for 1000 seconds with a sampling rate of 15 Hz. Only the last 2500 data points of the solution were used.

A.3 Coupled Rössler-Lorenz System

The coupled Lorenz-Rössler system is defined as

[TABLE]

where $b_{1}=8/3$, $b_{2}=0.2$, $c_{2}=5.7$, $k_{1}=0.1$, $k_{2}=0.1$, $k_{3}=0.1$, $\lambda=28$, $\sigma=10$, and $a=0.25$ for a periodic response and $a=0.51$ for a chaotic response. This system was simulated at a frequency of 50 Hz for 500 seconds with the last 300 seconds used.

A.4 Bi-Directional Coupled Rössler System

The Bi-directional Rössler system is defined as

[TABLE]

with $w_{1}=0.99$, $w_{2}=0.95$, and $k=0.05$. This was solved for 1000 seconds with a sampling rate of 10 Hz. Only the last 140 seconds of the solution were used.

A.5 Chua Circuit

Chua’s circuit is based on a non-linear circuit and is described as

[TABLE]

where $f(x)$ is based on a non-linear resistor model defined as

[TABLE]

The system parmeters were set to $\beta=27$, $\gamma=1$, $m_{0}=-3/7$, $m_{1}=3/7$, and $\alpha=10.8$ for a periodic response and $\alpha=12.8$ for a chaotic response. The system was simulated for 200 seconds at a rate of 50 Hz and the last 80 seconds were used.

A.6 Mackey-Glass Delayed Differential Equation

The Mackey-Glass Delayed Differential Equation is defined as

[TABLE]

with $\tau=2$, $\beta=2$, $\gamma=1$, and $n=9.65$. This was solved for 400 seconds with a sampling rate of 50 Hz. The solution was then downsampled to 5 Hz and the last 200 seconds were used.

A.7 Periodic Sinusoidal Function

The sinusoidal function is defined as

[TABLE]

This was solved for 40 seconds with a sampling rate of 50 Hz.

A.8 Quasiperiodic Function

This function is generated using two incommensurate periodic functions as

[TABLE]

This was sampled such that $t\in[0,100]$ at a rate of 50 Hz.

A.9 EEG Data

The EEG signal was taken from andrzejak et al. [3]. Specifically, the first 5000 data points from the EEG data of a healthy patient from set A (file Z-093) was used and the first 5000 data points of a patient experiencing a seizure from set E (file S-056) was used.

A.10 ECG Data

The Electrocardoagram (ECG) data was taken from SciPy’s misc.electrocardiogram data set. This ECG data was originally provided by the MIT-BIH Arrhythmia Database [30]. We used data points 3000 to 5500 during normal sinus rhythm and 8500 to 11000 during arrhythmia.

A.11 Logistic Map

The logistic map was generated as

[TABLE]

with $x_{0}=0.5$ and $r=3.95$. Equation 31 was solved for the first 500 data points.

A.12 Hénon Map

The Hénon map was solved as

[TABLE]

where $b=0.3$, $x_{0}=0.1$, $y_{0}=0.3$, and $a=1.4$. This system was solved for the first 500 data points of the x-solution.

A.13 Double Pendulum

The double pendulum is a staple benchtop experiment for investigated chaos in a mechanical system. A point-mass double pendulum’s equations of motion are defined as

[TABLE]

where the system parameters $g=9.81$ $\rm m/s^{2}$, $m_{1}=1$ kg, $m_{2}=1$ kg, $\ell_{1}=1$ m, and $\ell_{2}=1$ m. The system was solved for 200 seconds at a rate of 100 Hz and only the last 30 seconds were used with initial conditions $[\theta_{1},\theta_{2},\omega_{1},\omega_{2}]=[0,3{\rm\>rad},0,0]$. This system will have different dynamic states based on the initial conditions, which can vary from periodic, quasiperiodic, and chaotic.