Topological time-series analysis with delay-variant embedding

To qualitatively evaluate the characteristics of the time series, we apply topological data analysis Carlsson (2009), which is a computational method that can be used for characterizing the topological features of high-dimensional data. We construct a simplicial-complex model Edelsbrunner and Harer (2010) from the points in the embedded space and obtain the topological information as the number, position, and size of single- or multi-dimensional clusters and loops. We build a complex over a set of points if the pairwise distances between them are less than or equal to $2\varepsilon$, where $\varepsilon$ is a given non-negative scale parameter. Different values of $\varepsilon$ result in different complexes and different topological information. If $\varepsilon$ is considerably small, no connections are created; further, the resulting simplicial complex is not different when compared to the original points. As we gradually increase $\varepsilon$, connections appear between the points; however, if $\varepsilon$ becomes considerably large, all the points are connected with each other, and no useful information can be conveyed. If we increase $\varepsilon$ as $\varepsilon_{1}<\ldots<\varepsilon_{n}$, we obtain a sequence of embedded simplicial complexes that can be referred to as filtration (see Appendix A).

We use persistent homology theory Edelsbrunner et al. (2002); Zomorodian and Carlsson (2005) to study the topological features across filtration. A practical way to visualize the results of persistent homology is through multi-set points in the two-dimensional persistence diagram. In this diagram, each point $(b,d)$ represents an $l$-dimensional hole (i.e., the connected components are zero-dimensional, loops and tunnels are one-dimensional holes, and voids are two-dimensional holes) that appears at $\varepsilon=b$ (known as the birth scale) and disappears at $\varepsilon=d$ (known as the death scale) across the filtration (see Appendix A). This information encapsulates the topological features of the time series after embedding and provides valuable insights into the behavior of a dynamical system; for instance, the emergence of an oscillation in a time series can be attributed to the birth of a loop in the embedded space. Using these features, the qualitative properties of the time series can be captured in a robust and efficient manner.

For the observed time series $x(t)$, we use $F_{l,\tau}(x(t))$ to denote the two-dimensional persistence diagram calculated for $l$-dimensional holes from the embedded points with time delay $\tau$. We consider $\tau$ in a predefined set $\mathcal{T}=\{\tau_{1},\tau_{2},...,\tau_{K}\}$, where $\tau_{1}<\tau_{2}<\cdots<\tau_{K}$ ($K$ denotes the size of the set). The three-dimensional persistence diagram for time series $x(t)$ can be defined as

Figure 1 shows a schematic of the three-dimensional persistence diagram $\text{PD}^{(3)}_{1}(x(t))$ for loops and tunnels (one-dimensional holes), wherein the embedded points exhibit different shapes and topological features for different values of $\tau$. In the middle panel of Fig. 1, there are two big loops for $\tau_{1}$, whereas loops survive for short times with different distributions of birth and death scales (right panel) for $\tau_{2}$ and $\tau_{3}$. This information cannot be obtained using a single value of $\tau$ in the embedding.

Time-series data tend to be considerably noisy, which is a feature that can be considered to be either a measurement artifact or an inherent part of the dynamics themselves. Therefore, the persistence diagram should be stable with respect to the data being perturbed by noise. To evaluate the stability of the persistence diagram, we introduce the concept of bottleneck distance as a metric structure for comparing the persistence diagrams. Given two three-dimensional persistence diagrams $Dg_{1}^{(3)}$ and $Dg_{2}^{(3)}$, consider all matchings $\psi$ such that a point on one diagram can be matched either to a point on the other diagram or to its projection on the diagonal plane $\mathcal{W}^{(3)}=\{(b,b,\tau)\mid b,\tau\in\mathbb{R}\}$ (Fig. 2(a)). For each pair $(\bm{p},\bm{q})\in\psi$ for which $\bm{p}=\left(b_{1},d_{1},\tau_{1}\right)$ and $\bm{q}=\left(b_{2},d_{2},\tau_{2}\right)$, we define the relative infinity-norm distance between $\bm{p}$ and $\bm{q}$ as $d^{(\infty)}_{\xi}(\bm{p},\bm{q})=\max\left(|b_{1}-b_{2}|,|d_{1}-d_{2}|,\xi|\tau_{1}-\tau_{2}|\right)$, where $\xi$ is a positive rescaling coefficient introduced to adjust the scale difference between the point-wise distance and time. The bottleneck distance $d^{(3)}_{B,\xi}(Dg_{1}^{(3)},Dg_{2}^{(3)})$ can be defined as the infimum of the longest matched relative infinity-norm distance over all matchings $\psi$:

We show that three-dimensional persistence diagrams are stable with respect to the bottleneck distance under perturbation applied to the time series. Given two time series $x(t)$ and $y(t)$ of the same length, let $Dg_{(m)}^{(3)}(x)$ and $Dg_{(m)}^{(3)}(y)$ be their three-dimensional persistence diagrams respectively, as calculated for the embedding dimension $m$. Based on the stability properties of two-dimensional persistence diagrams Chazal et al. (2014), we can prove the following stability property of three-dimensional pesistence diagrams (see Appendix B):

for an arbitrary positive $\xi$. If we identify $y(t)$ as the perturbed data obtained by adding noise to $x(t)$, Eq. (3) shows that the upper limit of the bottleneck distance between $Dg_{(m)}^{(3)}(x)$ and $Dg_{(m)}^{(3)}(y)$ is governed by the magnitude of the noise. Thus, the inequality of Eq. (3) states that our three-dimensional persistence diagrams are robust with respect to the time-series data being perturbed by noise. Therefore, these diagrams can be used as discriminating features for characterizing the time series.

To use the persistence diagrams as features for performing the statistical learning tasks, such as classification, we must define a similarity measure, such as kernel mapping, for such diagrams Reininghaus et al. (2015); Kusano et al. (2016); Carrière et al. (2017). Because the points that are close to $\mathcal{W}^{(3)}$ can be considered to be insignificant topological features, they should not influence the computed value of this kernel. Given the positive bandwidth $\sigma$ and the positive rescaling parameter $\xi$, the kernel $k_{\sigma,\xi}$ between two three-dimensional persistence diagrams $Dg^{(3)}_{1}$ and $Dg^{(3)}_{2}$ can be defined as

where $\bm{\bar{q}}$ is a symmetric point of $\bm{q}$ with respect to $\mathcal{W}^{(3)}$ and $d^{2}_{\xi}(\bm{p},\bm{q})=|b_{1}-b_{2}|^{2}+|d_{1}-d_{2}|^{2}+\xi^{2}|\tau_{1}-\tau_{2}|^{2}$, $d^{2}_{\xi}(\bm{p},\bm{\bar{q}})=|b_{1}-d_{2}|^{2}+|d_{1}-b_{2}|^{2}+\xi^{2}|\tau_{1}-\tau_{2}|^{2}$, with $\bm{p}=\left(b_{1},d_{1},\tau_{1}\right)$ and $\bm{q}=\left(b_{2},d_{2},\tau_{2}\right)$. If $\bm{p}$ or $\bm{q}$ is near $\mathcal{W}^{(3)}$, i.e., $b_{1}\approx d_{1}$ or $b_{2}\approx d_{2}$, then $d^{2}_{\xi}(\bm{p},\bm{q})\approx d^{2}_{\xi}(\bm{p},\bm{\bar{q}})$ and the pair $(\bm{p},\bm{q})$ will have a minor influence on the value of the kernel (Fig. 2(b)). Based on Ref. Reininghaus et al. (2015), we can prove that $k_{\sigma,\xi}$ is a positive-definite kernel (see Appendix C). When we apply $k_{\sigma,\xi}$ to time-series classification, the parameters $\sigma$ and $\xi$ are observed to affect the classification performance and can be chosen either by cross-validation or in a heuristic manner Gretton et al. (2008) (see Appendix C). In our classification tasks, we use the normalized version of the kernel, which can be calculated as $k_{\sigma,\xi}(E,F)\leftarrow k_{\sigma,\xi}(E,F)/\sqrt{k_{\sigma,\xi}(E,E)k_{\sigma,\xi}(F,F)},$ where $E$ and $F$ denote two persistence diagrams. The source code that has been used to calculate the persistence diagrams and kernels can be found on GitHub Tran and Hasegawa .

Initially, the proposed method is applied to classify periodic and aperiodic time series using synthetic single-cell data. This is challenging because of the difficulty associated with discriminating between an oscillation containing noise and a mere noisy fluctuation Phillips et al. (2017). We generate synthetic mRNA and protein time-series data from a stochastic model of the Hes1 genetic oscillator Monk (2003); Galla (2009) exhibiting negative autoregulation with delay. We use the delayed version of the Gillespie algorithm Gillespie (1977); Anderson (2007) to generate data from 1,000 cells in both the oscillatory and non-oscillatory parameter regimes Galla (2009); Brett and Galla (2013). We measure the protein levels after every $\nu$ (= 64, 32, 16, 8) min for 4,096 min. We normalize the time series to have zero mean and unit variance; further, we add Gaussian white noise with variance $\sigma^{2}_{n}=0.1$ to assess the robustness of the method.

Figure 3 depicts the examples of the time series and persistence diagrams for the two regimes with an embedding dimension of $m=3$. Figures 3(a) and (d) show the time series generated by the measurements in the non-oscillatory and oscillatory regimes, respectively, after every $\nu=16$ min, and Figs. 3(b) and (e) denote their respective two-dimensional persistence diagrams computed from the single-delay embedding. In Figs. 3(b) and (e), $\tau$ is selected by the mutual-information method to maximize the statistical measure of nonlinear independence in the delay coordinates of the embedded points Fraser and Swinney (1986). Figures 3(c) and (f) depict the three-dimensional persistence diagrams obtained using Fig. 3(a) and (d), respectively, with $K=10$ time-delay values as $\tau_{1}=1,\tau_{2}=2,\ldots,\tau_{10}=10$. In these examples, it is difficult to distinguish between the two regimes using either the original time series or the two-dimensional diagrams, whereas there is a distinct pattern in the three-dimensional diagrams because all the topological variations are considered when $\tau$ changes. In the three-dimensional diagram of the non-oscillatory data, the points, which represent the loops in the embedded space, are extensively distributed along the birth and death scales, corresponding to the random fluctuation in the time series. In contrast, in the oscillatory data, points are densely distributed along the birth and death scales; this manifests the appearance of repeated loop patterns in the embedded space, thereby indicating periodic patterns in the time series.

Further, we assess the effectiveness of the proposed topological features by classifying the data obtained from the Hes1 model, which can be randomly split with equal probability to belong to the training and test sets. We use the support vector machine Bishop (2006) to perform classification in the kernel space. In the delay-variant method, we consider $K=10$ time-delay values as $\tau_{1}=1,\tau_{2}=2,\ldots,\tau_{10}=10$. In the single-delay method, we choose a single value of $\tau$ for each value from $\{\tau_{1},\ldots,\tau_{10}\}$. For each embedding dimension of the single-delay method, we use the maximum accuracy in the test dataset over all the values of $\tau$ for performing comparison with the delay-variant method. Figure 4 depicts the average classification accuracy obtained using 100 random splits at different embedding dimensions $m$ and different measurement intervals $\nu=8$ (Fig. 4(a)), $\nu=16$ (Fig. 4(c)), $\nu=32$ (Fig. 4(e)), and $\nu=64$ (Fig. 4(g)). When the same embedding dimension $m$ is being used, the total dimension for the delay-variant method ($mK$) is observed to be higher than that for the single-delay method ($m$). In the classification task, the high-dimensional features tend to achieve high classification accuracies. Therefore, to ensure that a fair comparison between the single-delay method and the delay-variant method can be performed, we define the effective embedding dimension as $M=mK$ and further compare the accuracy for different values of $M$ with $\nu=8$ (Fig. 4(b)), $\nu=16$ (Fig. 4(d)), $\nu=32$ (Fig. 4(f)), and $\nu=64$ (Fig. 4(h)). For $\nu=64$ (Fig. 4(h)), we do not consider embeddings with $M\geq 60$ because the length of the time series is 64.

When the comparison is performed over the embedding dimension $m$, Fig. 4 depicts that both the methods will achieve almost the same accuracy when the time series is long enough ($\nu=8,16$) at $m\geq 4$. In the single-delay method, $\tau$ should be selected by cross-validation for performing a fair comparison; however, even with trial and error to obtain the maximum accuracy in the test dataset in the single-delay method, our delay-variant method still performs better. When $\nu$ is increased, the time series is reduced, thereby increasing the difficulty of classification; regardless, the delay-variant method attains a higher accuracy than the single-delay method. For $\nu=32,64$, the delay-variant method outperforms the single-delay method over a wide range of $m$ ($4\text{--}10$) in terms of the accuracy. This observation indicates that the single-delay method degrades quickly with a shorter time series because it is highly sensitive to the choice of $\tau$. However, by considering a range of $\tau$ with delay embedding, this sensitivity can be reduced, thereby achieving stable performance.

While performing the comparison over the effective embedding dimension $M$, the accuracy of the delay-variant method does not change over a wide range of $M$ ($40\text{--}100$), indicating that this method is more reliable than the single-delay method (Fig. 4(b), (d), (f), and (h)). For a given effective dimension, the embedded points of single-delay method become highly folded and excessively sparse in a high-dimensional space because the embedding dimension that is used in the single-delay method is considerably higher than that used in the delay-variant method. Furthermore, selecting a considerably high embedding dimension increases the computational complexity of the embedding and causes the high impact of noise acting in a high proportion of elements in delay coordinates. Therefore, selecting a considerably high embedding dimension makes it disadvantageous to use topological features for characterizing the time series, thereby degrading the classification performance.

To demonstrate that the delay-variant method is more robust than the single-delay method when the time series is being perturbed by noise, we compare the two aforementioned methods based on the values of the noise variance $\sigma^{2}_{n}$ of 0.1, 0.3, 0.5, and 0.7. Figure 5 depicts the average classification accuracy versus the embedding dimension $m$ over 100 random splits at a measurement interval of $\nu=16$ min when the noise variance $\sigma^{2}_{n}=0.1$ (Fig. 5(a)), $\sigma^{2}_{n}=0.3$ (Fig. 5(b)), $\sigma^{2}_{n}=0.5$ (Fig. 5(c)), and $\sigma^{2}_{n}=0.7$ (Fig. 5(d)). As $\sigma^{2}_{n}$ is increased, the classification becomes more difficult because the time-series data become noisier. Even so, the delay-variant method still outperforms the single-delay method over $m$. Furthermore, the single-delay method degrades more quickly as $\sigma^{2}_{n}$ is increased, thereby indicating that our delay-variant method is more robust and stable in noisy conditions than the single-delay method.

Further, we evaluate the performance of the delay-variant method in classifying different heartbeat-signal patterns based on six real electrocardiogram (ECG) datasets in Chen et al. (2015), namely ECG200 (separate normal and myocardial infarction heartbeats; 200 time series of length 96), ECG5000 (separate five levels of congestive heart failure; 5,000 time series of length 140), ECGFiveDays (separate records from two different days for the same patient; 884 time series of length 136), TwoLeadECG (separate records from two different leads; 1,162 time series of length 82), and Non-InvasiveFetalECGThorax1 and Non-InvasiveFetalECGThorax2 (separate records from the left and right thorax with expert labeling of 42 classes of non-invasive fetal ECG; 3,765 time series of length 750 in each dataset). We use the delay-variant method to classify the Caenorhabditis elegans roundworms from EigenWorms dataset as either wild or mutant based on their movements Brown et al. (2013); Yemini et al. (2013). The movement trajectories are processed as 259 time series of length 900. Finally, we denote that our method can be used to diagnose whether a certain symptom can be observed in an automotive subsystem using the time-series data associated with the engine noise with the FordB dataset for 4,446 time series of length 500 Chen et al. (2015). Here, the training data in the FordB dataset were collected under typical operating conditions, whereas the test data were collected under noisy conditions. For these datasets, we employ the train-test split provided in Chen et al. (2015). For the Non-InvasiveFetalECGThorax1, Non-InvasiveFetalECGThorax2, EigenWorms, and FordB datasets, we downsampled the time series with sampling rates as 2, 2, 3, and 2, respectively, before computing the persistence diagrams.

In the delay-variant method, we use $K=10$ time-delay values of $\tau_{1}=1,\tau_{2}=2,\ldots,\tau_{10}=10$ with embedding dimensions of $m=2,3,\ldots,10$. In the single-delay method, we choose single $\tau$ for each value from $\{\tau_{1},\ldots,\tau_{10}\}$. For each embedding dimension of the single-delay method, we use the maximum accuracy in the test dataset over all values of $\tau$ to compare with the delay-variant method. In both the methods, we use a linear combination of normalized kernels for the zero-dimensional and one-dimensional holes. The combination weights and the embedding dimension $m$ are selected by cross-validation (see Appendix C). Further, we compare the delay-variant and single-delay methods along with some alternative standard approaches. Most of the previous research on time-series classification focused on finding appropriate similarity measures for the $1^{\text{st}}$-nearest neighbor (NN) classifier. For ensuring similarity in the time domain, we consider either the Euclidean (E) or dynamic time warping (D) distance. For ensuring similarity in the frequency domain and for similarity in the autocorrelation, we consider the power spectrum (PS) and the autocorrelation function (AC). We also examine the elastic ensemble (EE) Lines and Bagnall (2014) as a combination of nearest neighbor classifiers using multiple distance measures in the time domain. Finally, we perform the learned-shapelets (LS) method, which classifies time series by learning the representative shapelets (i.e., short discriminant time-series subsequences) Grabocka et al. (2014). Instead of implementing these standard algorithms, we use the results from Ref. Bagnall et al. (2016, ). The test results are presented in Table 1, wherein the best and second-best accuracy scores of each dataset are colored in dark pink and light pink, respectively. The delay-variant method shows better results than the single-delay method for all the datasets and outperforms all the other algorithms on an average. Further, the delay-variant method offers the best results for four of the eight datasets and the second-best results for the remaining four, suggesting that our method is an effective mechanism to classify the time series.

Subsequently, we investigate the types of features that can be effectively captured using the delay-variant method. Because the time scale of the extracted pattern is partially dependent on the time delay, the delay-variant method should be able to extract the patterns that comprise multiple different time scales. To better understand this ability, we analyze the synthetic noisy time series obtained from a frequency-modulated model. Given the original signal $s_{m}(t)=A_{m}\sin(2\pi f_{m}t)$ with a carrier signal $s_{c}(t)=A_{c}\sin(2\pi f_{c}t)$, we consider the modulated signal to be $s(t)=A_{c}\sin(2\pi f_{c}t+A_{m}\sin(2\pi f_{m}t))$, where $t\in[0,0.1]$, $A_{m}=10$, $A_{c}=1.0$, and $f_{m}=25$. Further, we consider the frequency $f_{c}$ to be $f_{c}=5f$ with $f=1\text{--}20$ to represent patterns with multiple time scales in the modulated signal. For each value of $f_{c}$, we generate 20 noisy discrete time series $z_{l}(n)=s(0.0002n)+w_{l}(n)$, where $l=1,2,\ldots,20$, $n=0,1,\ldots,500$, and $w_{l}(n)$ denote the Gaussian white noise with a variance of 0.1.

In the delay-variant method, we use $K=20$ time-delay values of $\tau_{1}=1,\tau_{2}=2,\ldots,\tau_{20}=20$. Figure 6 depicts a time series with $f_{c}=25$ and its three-dimensional persistence diagram for the loops or cycles at an embedding dimension of $m=6$. In the persistence diagram, the points near the diagonal plane $\mathcal{W}^{(3)}=\{(b,b,\tau)\mid b,\tau\in\mathbb{R}\}$ represent the noise in the time series, whereas the points that are away from $\mathcal{W}^{(3)}$ represent significant cycles in the embedded space. Therefore, the points in the diagram that has been plotted for different values of $\tau$ reveal the multiple-time-scale components in the time series corresponding to the limit-cycle oscillations of the dynamics.

An example of the principal component projection from the kernel space is depicted in Fig. 7 for (a) the delay-variant method with embedding dimension $m=6$, (b) the single-delay method with $m=6$ and $\tau=1,2,\ldots,20$, and (c) the single-delay method with $m=120$ (to compare with the delay-variant method at the same effective embedding dimension $M=120$). Because the length of the time series is 501, in the single-delay approach, $\tau$ can only take values from $\{1,2,3,4\}$ when $m=120$. Different colors represent the data generated using different values of $f_{c}$. Figure 7 depicts that the delay-variant method results in more distinct regions corresponding to different values of $f_{c}$ than that obtained using the single-delay method. Note that the points obtained using the single-delay method cannot be distinct, regardless of the angle. These results denote that the delay-variant method is superior to the single-delay method with respect to the ability to identify changes in the multiple-time-scale patterns in time series.