A Subsequence Approach to Topological Data Analysis for Irregularly-Spaced Time Series

Homology is an area of mathematics that looks for holes in a topological space, and persistent homology looks for holes in data. Topological holes are identified using ideas from algebraic topology and are characterized by different dimensional homology groups (Hatcher et al., 2002; Edelsbrunner and Harer, 2022). The zero-dimensional homology group (${H}_{0}$) contains connected components (clusters), the one-dimensional homology group (${H}_{1}$) contains loops, the two-dimensional homology group ($H_{2}$) contains voids like the interior of a balloon, and more generally, the k-dimensional homology group (${H}_{k}$) can be thought of as representing k-dimensional holes.

The data for which we use persistent homology in the proposed method are point clouds derived from a timeseries. If we only computed the homology of the point-cloud data, it would not be useful because the data points are not connected resulting in only $H_{0}$ features. Instead, an intermediate structure is defined with simplicial complexes. A zero-simplex is a point, a one-simplex is a segment, a two-simplex is a triangle, and, more generally, a $k$-simplex $\sigma=(v_{0},\cdots,v_{k})$ is a $k$-dimensional polytope of $k+1$ affinely independent points (i.e., zero-simplexes) $v_{0},\cdots,v_{k}$. A simplicial complex ${K}$ is a finite set of simplices such that for any $\sigma_{1},\sigma_{2}\in{K}$, $\sigma_{1}\cap\sigma_{2}$ is a face (i.e., the convex hull of any non-empty subset of points that define the simplex) of both simplices, or the empty set; and a face of any simplex $\sigma\in{K}$ is also a simplex in ${K}$. To compute the homology, simplicial complexes are built along a sequence of filtration values, which is discussed next.

A common approach to constructing simplicial complexes in TDA is the Vietoris-Rips (VR) complex (Vietoris, 1927; Edelsbrunner and Harer, 2022). A VR complex is constructed over a finite set of points $S=\{v_{1},v_{2},\cdots,v_{n}\}$ using a distance parameter $\delta\geq 0$. For any subset of k-points $\{v_{i_{1}},\cdots,v_{i_{k}}\}$, a $(k-1)$-dimensional simplex is formed when the pairwise Euclidean distance between all points is at most $\delta$. A collection of all such simplices forms the VR complex denoted as $VR(S,\delta)$. The composition of the VR complex progresses hierarchically as $\delta$ increases. This leads to the concept of filtration, which defines an inclusion relation between the simplicial complexes for a set of $\delta$ values. More formally, for an ordered sequence of $\delta$ values: $0<\delta_{1}<\delta_{2}<\cdots<\delta_{q}<\infty$, the VR complexes admits a nested structure, $VR(S,0)\subset VR(S,\delta_{1})\subset\cdots\subset VR(S,\delta_{q})\subset VR(S,\infty).$ Figure 1 provides an illustration of a VR complex. The black points (zero-simplices) in Figure 1(a) and Figure 1(b) denotes the data with cyan balls of diameter $\delta=0.8$ and $\delta=1.5$, respectively, along with the resulting one-simplices and two-simplices. The simplicial complex in Figure 1(a) is a subset of the simplicial complex displayed in Figure 1(b).

The inclusion relation between the VR complexes induces a map between the homology groups, ${H}_{k}(VR(S,0))\xrightarrow{}{H}_{k}(VR(S,\delta_{1}))\xrightarrow{}\cdots\xrightarrow{}{H}_{k}(VR(S,\infty)).$ The notion of persistent homology is developed through these homology maps by tracking the changes in the generators (i.e., the different dimensional holes) of these nested homology groups. The birth time and death time of homology group generators along this sequence encodes topological changes in the groups. For a homology group ${H}_{k}$, we denote the birth time and death time of the $j$-th homology group generator (i.e., the $j$-th hole) by $b_{j}$ and $d_{j}$, respectively. The persistence, or lifetime, of the generator is given by $d_{j}-b_{j}$, and a longer persistence is often considered to be topological signal while shorter persistence often represents topological noise. If we let $k_{j}$ be the homology group dimension of the $j$-th generator, and ${J}$ the index set of the generators of the homology groups, then the set $\text{Dgm}(S)=\{(b_{j},d_{j},k_{j}):\forall j\in J\}\cup\Delta$ fully characterizes the persistence (life span) of the homology generators, and is used to construct a graphical summary referred to as the persistence diagram; see example in Figure 1(c). The $\Delta$ represents points where the birth time is equal to the death time, and it included in the definition of persistence diagrams for technical reasons (Mileyko et al., 2011).

A persistence diagram is a summary statistic that characterizes the birth and death times of different dimensional holes detected in a dataset. These diagrams can be useful on their own for visualizing topological information of complex data, but are often used for inference or prediction tasks (e.g., Fasy et al., 2014; Bubenik et al., 2015; Adams et al., 2017; Robinson and Turner, 2017; Xu et al., 2019; Berry et al., 2020; Pun et al., 2022; Glenn et al., 2024). In this work, persistence diagrams are used to characterize the TDE from time series and to quantify the performance of the proposed and comparison methods for reconstructing the state space of irregularly-spaced time-series data.

For this work, the term uniformly-spaced time series is used to describe time series that have equally-spaced time intervals between successive observations, and is considered the “true” times series for purposes of evaluating the proposed method. The term irregularly-spaced time series refers to observations with unequally-spaced time interval between successive observations. To characterize how the two forms of time series are related, it is assumed throughout this work that the irregularly-spaced time series is a subset of the uniformly-spaced time series. Hence the irregularly-spaced time series always have fewer time measurements than the corresponding uniformly-spaced time series.

The concept of TDE is also referred to in the literature as a delayed coordinate embedding, a sliding window embedding, or simply a Takens embedding. For this work, only the term TDE is used. For any irregularly-spaced time series, it is assumed there is a “true” underlying uniformly-spaced time series. The use of “TDE” exclusively refers to an embedding constructed from this true underlying uniformly-spaced time series. After applying the proposed subsequence method, the resulting embedding is referred as the subsequence embedding (SSE). In instances where an exposition applies to both the TDE and SSE, the term embedding map is used as a collective reference to the two concepts.

For a given time interval, it is assumed that missing or unobserved values occurs with a given probability. That is, for a given time point in a time interval, a measurement is not observed at that point with some probability. Such probabilistic mechanism governing the observations of time series values is not uncommon in the literature (e.g., Dunsmuir and Robinson, 1981). This probability can be fixed for all time points or it can vary for each time point. This characterization is referred to as the missingness structure of the time series in context.

For this work, the discussion on TDEs is restricted to univariate time series. Assume this univariate time series is generated by a system with a state vector $\mathbf{s}(t)$ on a manifold which is a subset of some $N$-dimensional space $\mathbb{R}^{N}$. The state vector $\mathbf{s}(t)$ is not directly observable, however some measurement of it, denoted $x(t)=h(\mathbf{s}(t))$, is observed through the measurement function $h(\cdot)$. The measurement function $h(\cdot)$ can be thought of as rule that transforms the high-dimensional state vector into the observed univariate time series $x(t)$. For instance, in astronomy, $\mathbf{s}(t)$ might include variables such as positions, velocities, and luminosities of various celestial bodies, such as exoplanets, stars, or galaxies. However, the measurement function is specifically designed to extract a single scalar value from this vector. The specific form of $h(\cdot)$ is influenced by many considerations, for example, the limitation of observational tools. For a star, the measurement function could be designed to extract a key observable from the state vector, such as its luminosity. Thus the measurement $h(\mathbf{s}(t))$ reflects the observed brightness of the star at any given time $t$.

The scalar value $x(t)$ is the observed time series measurement. Define the function $F:\mathbb{R}^{N}\xrightarrow{}\mathbb{R}^{M+1}$ as the embedding map with the form $F(\mathbf{s}(t))=\left[x(t),x(t+\tau),\cdots,x(t+M\tau)\right]$. If the measurement function $h(\cdot)$ is noise-free, and the embedding dimension $M+1$ is chosen to be more than twice the dimension of the attractor (i.e., the $N$-dimensional region toward which the system evolves) of the system’s state space, Takens’s theorem guarantees that the embedding map $F(\mathbf{s}(t))$ has a one-to-one correspondence between the original state space of the system (from which the time series is derived) and the reconstructed state space formed by $F(\mathbf{s}(t))$ (Takens, 2006). This ensures that the dynamics of the system can be studied in the reconstructed space as if it were being studied in the original space. Figure 2 demonstrates this reconstruction process by mapping the scalar time series to the TDE matrix $F$ to reconstruct the state space.

The choice of embedding dimension $M+1$ and step size $\tau$ is a subject of considerable research in the literature (e.g., Cao, 1997; Kim et al., 1999). In this work, the embedding dimension is chosen manually. However, the method of false nearest neighbors is one common method for determining this dimension, which identifies points in a low-dimensional space that appear to be near each other but are not actually neighbors when the data is viewed in a high-dimensional space. By systematically increasing the embedding dimension, and evaluating the percentage of false nearest neighbors, the dimension can be set where this percentage drops significantly, indicating a suitable dimension. More details about this and other procedures for determining $M$ and $\tau$ can be found in Cao, (1997); Kim et al., (1999). A large value of $M$ is often preferred as it enables the embedding to capture more details inherent in the time series. If $M$ is too large, there may be an insufficient number of points in the embedding space. Furthermore, if $M\tau$ is too small due to a small $\tau$, relatively fewer points fall in each embedding window. This results in points repeatedly appearing in windows, which can lead to redundant information. If $M\tau$ is too large due to a large value of $\tau$, the reconstructed state space can be distorted because relevant periodic behavior of the time series may not be captured (Casdagli et al., 1991). Hence the choice of $\tau$ and $M$ is such that $M\tau$ is not too large or too small, but is application-dependent and requires empirical testing.

The TDE construction in the previous section assumes the observed time series is uniformly-spaced, but a time series is often irregularly-spaced in real data. We propose a method to extract uniformly-spaced subsequences from the observed irregularly-spaced time series and prove its topology-preserving properties, along with consistency and convergence results.

The set of persistence diagrams from the true state space (TDE) and the reconstructed state space (SSE) can be endowed with a distance measure, such as the popular bottleneck distance. The bottleneck distance gives the minimal bijection between any two diagrams as a measure of the distance between the diagrams. Let $\mathbf{F}^{1}$ and $\mathbf{F}^{2}$ be two finite compact subsets of $\mathbb{R}^{M+1}$, with $\text{Dgm}(\mathbf{F}^{1})$ and $\text{Dgm}(\mathbf{F}^{2})$ as their corresponding persistence diagrams. The bottleneck distance, $\text{d}_{\text{B}}$, between the two persistence diagrams is defined as

where the infimum is taken over all bijections $\gamma:\text{Dgm}(\mathbf{F}^{1})\xrightarrow{}\text{Dgm}(\mathbf{F}^{2})$. For the metric space $(\mathbb{R}^{M+1},\lVert\cdot\rVert_{2})$, where $\lVert\cdot\rVert_{2}$ denotes the $l^{2}$-norm, define the distance function as: $d(F,\mathbf{F}^{1})=\inf_{F^{1}\in\mathbf{F}^{1}}\lVert F-F^{1}\rVert_{2},$ $\forall F\in\mathbb{R}^{M+1}.$ Then the Hausdorff distance, $d_{H}$, is given by

The Hausdorff distance quantifies the level to which two metric spaces represent the same object. A modification of the Hausdorff distance that is blind to isometric transformations is the Gromov-Hausdorff distance. Let $\nu_{1}:\mathbf{F}^{1}\hookrightarrow\mathbb{R}^{M+1}$ and $\nu_{2}:\mathbf{F}^{2}\hookrightarrow\mathbb{R}^{M+1}$ be isometric (distance preserving) embeddings into the metric space $(\mathbb{R}^{M+1},\lVert\cdot\rVert_{2})$. Then the Gromov-Hausdorff distance between $\mathbf{F}^{1}$ and $\mathbf{F}^{2}$ is given by: $d_{GH}(\mathbf{F}^{1},\mathbf{F}^{2})=\inf_{\nu_{1},\nu_{2}}d_{H}(\nu_{1}(\mathbf{F}^{1}),\nu_{2}(\mathbf{F}^{2}))$. A fundamental result on persistence diagrams is that they are stable summaries in many settings (i.e., a small change in a point cloud result in a small change in the corresponding persistence diagram) (Chazal and Michel, 2021). The stability relation of the persistence diagrams is stated as

In what follows, these stability results are established for the proposed SSE method and a denoising procedure to reduce the noise present in the observed time series.

Time-series data is typically observed with noise. The level of noise in the reconstructed space influences the presence and the persistence of the topological features. A Fourier denoising procedure is proposed to filter out noise in the observed time series.

Let $\mathbf{x}=\begin{bmatrix}x(t_{1}),&x(t_{2}),&\cdots&,x(t_{n})\end{bmatrix}^{\top}$ be an observed time series vector. The first step in the denoising procedure is to transform this observed signal to the frequency domain. The discrete Fourier transform (DFT) of $x(t_{k})$, denoted as $\tilde{x}(t_{k})$ is given by

where $j$ is the imaginary unit ($j^{2}=-1$), $\phi_{kr}=e^{-j2\pi w_{r}f_{k}}$, $0\leq w_{r}\leq 1$ are sample points, $0\leq f_{k}\leq n$ are frequencies, and $\tilde{x}(t_{k})$ is the $k$-th sample of the power spectrum at $f_{k}$.

To filter out noise, the power spectral density is computed for each observation $\tilde{x}(t_{k})$. Then a threshold is chosen, and any $\tilde{x}(t_{k})$ with power spectral density less than the threshold is set to zero. In selecting the threshold, the goal is to choose a value that does not smooth out the peaks in the true signal. The derivation that follows assumes the selected threshold preserves the peaks in the true signal. To simplify notations, the thresholded observations are also denoted as $\tilde{x}(t_{k})$. The thresholded $\tilde{x}(t_{k})$ are transformed back to the time domain to get the noise-reduced signal. which typically involves multiplying $\tilde{x}(t_{k})$ by the inverse of a Fourier transform matrix.

Let $\tilde{\mathbf{x}}$ be the DFT of the time series vector with corresponding Fourier basis $\boldsymbol{\phi}_{k}$, where

The forward transform in Equation (6) can be vectorized:

The backward transform can then be determined by inverting the matrix $\boldsymbol{\Phi}$. However, due to the non-uniformity in the spacing of the time series $\mathbf{x}$, the columns of $\boldsymbol{\Phi}$ are not orthogonal, and it is not directly invertible, so the pseudo-inverse is used instead. The backward transform then has the form: $\mathbf{x}=\frac{1}{n}(\boldsymbol{\Phi}^{H}\boldsymbol{\Phi})^{\dagger}\boldsymbol{\Phi}^{H}\tilde{\mathbf{x}}$, where $\mathbf{A}^{H}$ and $\mathbf{A}^{\dagger}$ denote the complex conjugate transpose and the Moore-Penrose inverse of the matrix $\mathbf{A}$, respectively. The matrix $(\boldsymbol{\Phi}^{H}\boldsymbol{\Phi})^{\dagger}\boldsymbol{\Phi}^{H}$ projects the frequency vector $\mathbf{\tilde{x}}$ onto the column space of $\boldsymbol{\Phi}$. Let $\Pi_{\boldsymbol{\Phi}}(\mathbf{x})$ denote this projection operation. The modulus of off-diagonal elements of $\boldsymbol{\Phi}^{H}\boldsymbol{\Phi}$ is bounded as: $\left|\sum_{k=1}^{n}e^{j2\pi(w_{l_{1}}-w_{l_{2}})f_{k}}\right|\leq\sum_{k=1}^{n}\left|e^{j2\pi(w_{l_{1}}-w_{l_{2}})f_{k}}\right|=n$, where the equality follows from the definition of the complex modulus $|{z}|=\sqrt{{z}\bar{z}}$, with $\bar{z}$ as the conjugate of the complex value $z$. The matrix $\boldsymbol{\Phi}\boldsymbol{\Phi}^{H}$ has the structure

for indicator function $\mathbbm{1}(l_{1}\leq l_{2})$ , and $(\boldsymbol{\Phi}\boldsymbol{\Phi})^{H}_{l_{1},l_{2}}$ denotes the value in the $l_{1}$-th row and $l_{2}$-th column of $\boldsymbol{\Phi}\boldsymbol{\Phi}^{H}$. Then $\boldsymbol{\Phi}\boldsymbol{\Phi}^{H}$ is Toeplitz when the frequency components are uniformly-spaced such that $f_{k}=k$, and $\boldsymbol{\Phi}\boldsymbol{\Phi}^{H}$ is fully specified by its first row elements (HuoLiu and YuanTang, 1998).

Using these results, the following proposition asserts that the DFT preserves topological features and is stable with respect to the bottleneck distance. In particular, the bottleneck distance between the persistence diagrams of the embeddings of the noise-free and smoothed (i.e., noise-reduced) time series is bounded above by the embeddings of the observed noisy and noise-free time series.

In this section, an error bound and stability results for the SSE approximation of the TDE are established. To simplify the notation, the embedding matrices are represented as sets where the elements of the set are the row vectors of the corresponding TDE or SSE. Also, for an embedding from a single uniformly-spaced time series, the step-size is assumed to be $\tau$. When constructing from a set of $P$ subsequences, a step-size of $\tau_{p}$ is assumed for the $p$-th subsequence, where $1\leq p\leq P$.

Because the SSE can have fewer elements than the TDE, the following lemma addresses how to expand the SSE without affecting its topology by repeating already existing points in the SSE, so distances can be computed between the TDE and the expanded SSE.

Lemma 4.2 asserts that a particular set of points can be added to an embedding matrix without changing the SSE’s persistence diagram. This is used to establish a bound on the SSE as an approximation to the TDE in the following proposition. In Lemma 4.2, when $k=(n-M\tau)-\sum_{p=1}^{P}(n_{p}-M\tau_{p})$, the row dimension of $\widehat{\mathbf{F}}^{2}$ is the same as that of $\mathbf{F}^{1}$. In such instances, when the interest is in a row-wise comparison of $\widehat{\mathbf{F}}^{2}$ and $\mathbf{F}^{1}$ the subsequence indexing in the time variable for any $\widehat{F}^{2}(\mathbf{s}(t_{p,i}))\in\widehat{\mathbf{F}}^{2}$ is ignored and a row is simply written as $\widehat{F}^{2}(\mathbf{s}(t_{i}))$, where $1\leq i\leq n-M\tau$.