Visualizing High Dimensional Dynamical Processes

TL;DR

This paper introduces DIG, a visualization method for multivariate time series that uses information geometry derived from diffusion operators to reveal complex data structures, demonstrated on EEG sleep stage data.

Contribution

The paper presents a novel visualization technique, DIG, which employs a new group of diffusion-based distances to better reveal data structure in high-dimensional dynamical processes.

Findings

DIG effectively visualizes sleep stages in EEG data.

New diffusion-based distances reveal structures not seen with traditional methods.

Method enhances understanding of high-dimensional dynamical systems.

Abstract

Manifold learning techniques for dynamical systems and time series have shown their utility for a broad spectrum of applications in recent years. While these methods are effective at learning a low-dimensional representation, they are often insufficient for visualizing the global and local structure of the data. In this paper, we present DIG (Dynamical Information Geometry), a visualization method for multivariate time series data that extracts an information geometry from a diffusion framework. Specifically, we implement a novel group of distances in the context of diffusion operators, which may be useful to reveal structure in the data that may not be accessible by the commonly used diffusion distances. Finally, we present a case study applying our visualization tool to EEG data to visualize sleep stages.

Equations

$$K({\mathbf{x}}_i,{\mathbf{x}}_j)=\exp \left(-\frac{||{\mathbf{x}}_i-{\mathbf{x}}_j|{|}^2}{2{\sigma }^2}\right),$$

$$D_t({\mathbf{x}}_i,{\mathbf{x}}_j{)}^2=\sum _{m=1}^N\frac{([P^t{]}_{mi}-[P^t{]}_{mj}{)}^2}{{\varphi }_0(m)},$$

$$K_{k,\alpha }({\mathbf{x}}_i,{\mathbf{x}}_j)=\frac{1}{2}\exp {\left(-\frac{||{\mathbf{x}}_i-{\mathbf{x}}_j|{|}^2}{{\sigma }_k({\mathbf{x}}_i)}\right)}^{\alpha }+\frac{1}{2}\exp {\left(-\frac{||{\mathbf{x}}_i-{\mathbf{x}}_j|{|}^2}{{\sigma }_k({\mathbf{x}}_j)}\right)}^{\alpha }.$$

$${\mathbf{z}}_t$$

$$d{\theta }_t^i$$

$$d^2({\mathbf{z}}_t,{\mathbf{z}}_s)=(\mathbb{E}({\mathbf{h}}_t)-\mathbb{E}({\mathbf{h}}_s){)}^T({\mathbf{C}}_t+{\mathbf{C}}_s{)}^{-1}(\mathbb{E}({\mathbf{h}}_t)-\mathbb{E}({\mathbf{h}}_s)),$$

$$\Vert {\mathbf{\theta }}_t-{\mathbf{\theta }}_s{\Vert }^2\approx d^2({\mathbf{z}}_t,{\mathbf{z}}_s).$$

$$d^2({\mathbf{z}}_t,{\mathbf{z}}_s)=\frac{1}{2}\sum _{i=1}^N\text{ln}({\lambda }_i),$$

$$\displaystyle{\scriptstyle D^{\gamma}_{t}(\bm{z}_{i},\bm{z}_{j})^{2}}=\left\{\begin{array}[]{ll}\sum_{m=1}^{N}\frac{(\text{log}[P^{t}]_{mi}-\text{log}[P^{t}]_{mj})^{2}}{\phi_{0}(m)},&\scriptstyle{\gamma=1}\\ \sum_{m=1}^{N}\frac{([P^{t}]_{mi}-[P^{t}]_{mj})^{2}}{\phi_{0}(m)},&\scriptstyle{\gamma=-1}\\ \sum_{m=1}^{N}\frac{2(([P^{t}]_{mi})^{\frac{1-\gamma}{2}}-([P^{t}]_{mj})^{\frac{1-\gamma}{2}})^{2}}{(1-\gamma)\phi_{0}(m)},&\scriptstyle{-1<\gamma<1}.\\ \end{array}\right.$$

$$D({\mathbf{z}}_i,{\mathbf{z}}_j)=2{\text{cos}}^{-1}\left(\sum _{m=1}^N\sqrt{[P^t{]}_{mi}[P^t{]}_{mj}}\right).$$

$$\mathtt{Stress}({\hat{z}}_1,\dots ,{\hat{z}}_N)=\sqrt{\frac{{\sum }_{i,j}{\left(D({\mathbf{z}}_i,{\mathbf{z}}_j)-\Vert {\hat{x}}_i-{\hat{x}}_j\Vert \right)}^2}{{\sum }_{i,j}{\left(D({\mathbf{z}}_i,{\mathbf{z}}_j)\right)}^2}},$$

Full text

Index Terms—  Visualization, dynamical processes, EIG, PHATE, diffusion maps

1 Introduction

Manifold learning techniques have become of great interest when studying high dimensional data. The underlying idea behind these methods, is that high dimensional data often encapsulates redundant information. In these cases, the data have an extrinsic dimensionality that is artificially high, while its intrinsic structure is well-modeled as a low-dimensional manifold plus noise. Following the same line of reasoning, dynamical systems and time series can be regarded as processes governed by few underlying parameters, confined in a low-dimensional manifold [1, 2].

In particular, electroencephalographic (EEG) measures can be contemplated in this analytical framework. These measures are taken from different parts of the brain, resulting in a multivariate time series in a high dimensional space. It is known that these time series are highly correlated with each other. Therefore, it can be assumed that there is a low-dimensional representation of the intrinsic dynamics of the brain that can explain a broad spectrum of physical and psychological phenomena such as sleep stages. Additionally, it can be very useful for researchers to achieve meaningful visual representations of this phenomena in two or three dimensions. Such visualizations can be used to better understand the overall shape and finer patterns within the data.

In this paper we present DIG (dynamical information geometry), a dimensionality reduction tool that is designed for visualizing the inherent low-dimensional structure present in high-dimensional dynamical processes. DIG is built upon a diffusion framework adapted to dynamical processes, followed by an embedding of a novel group of information distances applied to the diffusion operator. The resulting embedding is noise resilient and presents a faithful visualization of the true structure at both local and global scales with respect to time and the overall structure of the data. We demonstrate our DIG on high-dimensional EEG data.

This paper is organized as follows. In Section 2 we present a brief background of diffusion methods. In Section 3, we present the steps of DIG, with a focus on the extension of the diffusion process to dynamical systems and the embedding of information distances. Then, in Section 4, we show and discuss some practical results for EEG data. And finally in Section 5, we conclude our work and propose some extensions.

2 Background

2.1 Related Work

Many dimensionality reduction methods exist, some of which have been used for visualization [3, 4, 5, 6, 7, 8, 9]. Principal components analysis (PCA) [8] and t-distributed stochastic neighborhood embedding (t-SNE) [3] are two of the most commonly used methods for visualization. However, these and other methods are inadequate in many applications. First, these methods tend to favor one aspect of the data at the expense of the other. For example, when used for visualization, PCA typically shows the large scale global structure of the data while neglecting the finer, local structure. In contrast, t-SNE is explicitly designed to focus on the local structure and often distorts the global structure, potentially leading to misinterpretations [10]. Second, PCA and t-SNE fail to explicitly denoise the data for visualization. Thus in noisy settings, the true structure of the data can be obscured. In addition, none of these methods are designed to exploit the structure present in dynamical systems.

Diffusion maps (DM) is a popular nonlinear dimensionality reduction technique that effectively denoises the data while capturing both local and global structure [11]. However, DM typically encodes the information in higher dimensions and is thus not optimized for visualization. A more recent visualization method called PHATE was introduced in [12] to exploit the power of DM in denoising and capturing the structure of data while presenting the learned structure in low-dimensions by preserving an information distance between the diffusion probabilities.

DM has been extended to dynamical systems previously [13, 14, 15, 16]. In particular, Talmon and Coifman [14, 15] introduced an approach called empirical intrinsic geometry (EIG) that builds a diffusion geometry using a noise resilient distance. The resulting embedding learned from the geometry is thus noise-free and captures the true structure of the underlying process. However, EIG and other extensions of DM to dynamical systems are still not optimized for visualization as the learned structure of the data is encoded in higher dimensions. In this work, we introduce a new visualization method DIG that is well-suited for visualizing high-dimensional dynamical processes by preserving an information distance between the diffusion probabilities constructed from a noise resilient distance. This results in a visualization method that represents the true structure of the underlying dynamical process.

EEG signals have been embedded in low dimensional representations for detecting emotional states [17], pre-seizure states [2, 18] and sleep dynamics [16]. In the latter DM is implemented by building the affinity matrix using both the cross-spectrum distance and the covariance matrix distance as similarity measures between multivariate time series. EIG has also been applied to data including both respiratory and EEG signals [19].

2.2 Preliminaries

Here we provide background on DM [11] and PHATE [12]. DM learns the geometry of the data by first constructing a graph based on local similarities. The graph is typically constructed by applying a Gaussian kernel to Euclidean distances:

[TABLE]

where $\bm{x}_{1},\dots,\bm{x}_{N}$ are the data and $\sigma$ is a fixed kernel scale or bandwidth that controls the locality scale of the graph.

The resulting kernel or affinity matrix $K$ encodes the local structure of the data. The process of diffusion is then used to learn the global structure and denoise the data. The first step is to transform the affinity matrix $K$ into a probability transition matrix $P$ (also known as the diffusion operator), by dividing each row by the sum of its entries. The $i,j$th entry of $P$ indicates the probability of transitioning from $\bm{x}_{i}$ to $\bm{x}_{j}$ in a single-step random walk, where the probabilities are calculated based on the relative affinity between points. Hence, if the affinity between two points is high, then the transition probability is high. The global structure of the dynamical process is then learned by performing a $t$-step random walk for integers $t\geq 1$. The transition probabilities of these random walks are given by the $t$-th power of the diffusion operator, (i.e., $P^{t}$), and rows of this matrix serve as $t$-step diffusion representations of data points.

In DM, the information encoded in the diffused operator $P^{t}$ is typically embedded in lower dimensions via eigendecomposition. The Euclidean distances between the embedded coordinates are equivalent to the following scaled distance between entries in the diffused operator:

[TABLE]

where $\phi_{0}$ is the stationary distribution of the corresponding Markov chain, or equivalently the first left eigenvector of $P$.

DM, as described above, has several weaknesses. First, in many applications the data are not sampled uniformly. In these cases, a fixed bandwidth for all points with the Gausssian kernel in (1) may not accurately capture the local data geometry in all settings. For example, a bandwidth tuned for a densely sampled region will be inappropriate for sparsely sampled regions and vice versa. PHATE counters this by replacing the Gaussian kernel with fixed bandwidth with the $\alpha$-decay kernel with adaptive bandwidth [12], eq. (3). The adaptive bandwidth enables the affinities to scale with the local density while the $\alpha$-decay kernel corrects inaccuracies that may be introduced by the adaptive bandwidth in sparse sampled regions:

[TABLE]

The value of $\sigma_{k}(\bm{x}_{i})$ is the distance from $\bm{x}_{i}$ to its $k$-nearest neighbor. Since this value changes for different observations, it may result in a non-symmetric kernel. Taking the average, as shown in eq. (3), mitigates this issue.

Second, choosing an appropriate time scale $t$ is a difficult problem, as this parameter controls the resolution of the captured diffusion representation. In PHATE, the von Neumann entropy (VNE) is used to automatically tune this scale. The VNE (also known as spectral entropy) for each $t$ is the entropy of eigenvalues of the diffused operator $P^{t}$ (normalized to form a probability distribution), and provides a soft proxy for the number of significant eigenvalues of $P^{t}$. As $t\rightarrow\infty$, the VNE converges to zero, since the diffusion process is constructed to have a unique stationary distribution. The rate of decay of the VNE as $t$ increases is thus used to determine the appropriate value of $t$ at the transition between rapid decay, which is interpreted in [12] as corresponding to the elimination of noise, and slow decay (interpreted there as losing meaningful information).

Third, DM is not well-suited for visualization as the eigendecomposition does not explicitly force as much information as possible into low dimensions. In fact, when the data have a branching structure, DM tends to place different branches in different dimensions [20]. Additionally, attempts to directly embed the diffusion distances into low dimensions for visualization using multidimensional scaling (MDS) [7] can lead to unstable and inaccurate embeddings [12]. PHATE counters this by constructing a potential distance (i.e., comparing energy potentials) from diffused probabilities and then directly embedding potential distances in two or three dimensions using MDS.

3 The DIG Algorithm

In this section, we extend principles of DM and PHATE to dynamical systems to derive DIG. In this context, we present a family of information distances and derive some of their properties.

3.1 Diffusion with Dynamical Systems

In the context of dynamical systems it is needed to learn the local structure by constructing a matrix that encodes the local distances between data points. These local distances can be taken as an input for (3) to build a diffusion operator, from which information is extracted for visualization. To do this, we build upon the EIG framework [15, 14] which uses a state-space formalism (4)-(5):

[TABLE]

The multivariate time series $\bm{z}_{t}$ represents the observed time series data while $\bm{\theta}_{t}$ represents the hidden (unobserved) states that drive the process. $\bm{z}_{t}$ can be viewed as a corrupted version of a clean process $\bm{y}_{t}$ that is driven by the hidden states, where the corruption $\bm{\xi}_{t}$ is a stationary process independent of $\bm{y}_{t}$. In general, we can view $\bm{y}_{t}$ as being drawn from a conditional pdf $p(\bm{y}|\bm{\theta})$. In the stochastic process (5), the unknown drift functions $a^{i}$ are independent from $\theta^{j}$, $j\neq i$. Therefore, we assume local independence between $\theta^{i}_{t}$ and $\theta^{j}_{t}$, $\forall i\neq j$. The variables $w^{i}_{t}$ are Brownian motions.

It can be shown that the pdf $p(\bm{z}|\bm{\theta})$ is a linear transformation of $p(\bm{y}|\bm{\theta})$ [15, 14]. Since the pdfs are unknown, we use histograms as their estimators. Each histogram $\bm{h_{t}}=(h_{t}^{1},\ldots,h_{t}^{Nb})$ has $Nb$ bins, and is built with the observations within a time window of length $L_{1}$, centered at $\bm{z_{t}}$. The expected value of the histograms, e.g. $\mathbb{E}(h_{t}^{j})$, is a linear transformation of $p(\bm{z}|\bm{\theta})$. Since the Mahalanobis distance is invariant under linear transformations, it can be deduced that the distance (6) is noise resilient [14]:

[TABLE]

where $C_{t}$ and $C_{s}$ are the covariance matrices in the histograms space, in a time window of length $L_{2}$, centered at $\bm{h}_{t}$ and $\bm{h}_{s}$, respectively. Also, it can be proved under certain assumptions that $d(\bm{z}_{t},\bm{z}_{s})$ is a good approximation of the distance between the underlying state variables [14]:

[TABLE]

To learn the global relationships from the local information encoded by the distances in eq. (6), we use the diffusion process. We apply the $\alpha$-decay kernel in eq. (3) to the Mahalanobis distances in eq. (6), where the $k$-nearest neighbor distances are also calculated with this distance. We use the $\alpha$-decay kernel to account for potentially different regions of density, as described in Section 2.2. The diffusion operator $P$ is then constructed by row-normalizing the resulting kernel matrix as before.

For comparison purposes, we also make use of an alternative distance to (6). Assuming that the data within time windows of length $L_{1}$ centered at $z_{t}$ follows a multivariate Gaussian distribution $\mathcal{N}(\bm{\mu},\Sigma_{t})$, we can compute the geodesic distance between different time windows of data using the Fisher information as the Riemannian metric [21] as follows:

[TABLE]

where $\lambda_{i}$ are the roots of, $|\Sigma_{t}-\lambda\Sigma_{s}|=0$. The diffusion operator can then be obtained using this distance as input to the $\alpha$-decay kernel.

3.2 Embedding Information Distances

Typically, information is extracted from the diffusion operator $P$ by either eigendecomposition or by embedding the diffusion distances. However, as mentioned previously, the former typically fails to provide a low-dimensional representation that is sufficient for visualization while the latter can result in unstable embeddings in some cases [12]. To overcome this, DIG extracts the information from the diffusion operator by embedding an information distance instead. We focus on a broad family of information distances that are parametrized by $\gamma$:

[TABLE]

The parameter $\gamma$ controls the level of influence of the lower differences among probabilities in the overall distance. For example, the standard diffusion distances ($\gamma=-1$) are highly influenced by the highest absolute differences among probabilities. In contrast, the potential distances ($\gamma=1$), which were used in PHATE, account for the relative differences between them. Thus the standard diffusion distances and the potential distances can be viewed as two extremes of a general class of distances over the diffusion geometry.

It can be shown that for $\gamma\in[-1,1]$, the distance $D_{t}^{\gamma}$ forms an M-divergence [22, 23]. Furthermore, when $\gamma=0$, $D_{t}^{\gamma}$ becomes proportional to the Hellinger distance, which is an $f$-divergence [24, 25]. Further, $f$-divergences are directly related to the Fisher information and thus are well-suited for building an information geometry [26]. Therefore, $f$-divergences may be desirable for embedding the diffused probabilities.

We also consider another information distance based on $f$-divergences that has not been applied to diffusion operators as far as we know. Since the rows of the diffusion matrix $P$ can be interpreted as multinomial distributions, we can compute the geodesic distance between them using the Fisher information as the Reimmanian metric [27]. This can be seen as an extension of the KL divergence or the Hellinger distance, both $f$-divergences, for distributions far apart from each other. Such distances are as follows [28, 21]:

[TABLE]

After the information distances have been obtained, DIG applies metric MDS to the information distances to obtain a low-dimensional representation. Given an information distance $D$, metric MDS minimizes the following stress function:

[TABLE]

where the $\hat{z}_{i}$ are the $m$-dimensional embedded coordinates. For visualization, $m$ is chosen to be 2 or 3. See Algorithm 1 for pseudocode summarizing the described steps of DIG.

4 Experimental results

4.1 Data

We now present a real-world data application using EEG data provided by [29, 30]. The original data is sampled at 512Hz and labeled for every 30 second interval, within six sleep categories according to R&K rules (REM, Awake, S-1, S-2, S-3, S-4). Due to the lack of observations in some stages, we group S-1 with S-2, and S-3 with S-4. We band-filtered the data between 8-40 Hz, and down-sampled it to 128Hz.

4.2 Experimental Setup

The tuning and influence of the parameters $\alpha$, $t$ and $k$ for nondynamical data have been covered in [12]. For the EEG data, we found that the visualizations are highly robust for a wide range of configurations. Preliminary experiments showed that $k=5$, $\alpha=10$, and $t=10$ give meaningful results. Therefore, we used these parameters for all experiments shown here unless otherwise stated. To compute the distances (6) and (8), we need to choose $L_{1}$, the window size. Its selection is driven by the way the data is presented. In our case, we simply took $L_{1}=3840$, the number of observations in the 30s span. We also selected $L_{2}=10$, and the number of bins for the histograms $Nb=20$. We do not focus on these parameters since their impact have been already studied in the cited literature.

Conversely, as far as our knowledge goes, there has not been any previous study addressing the properties of the information distances mentioned in the previous section. Thus, we focus on how these distances may affect the learning process.

For this purpose, we wish to measure the relative distortion of the embeddings both locally and globally, when setting different values of $\gamma$, as well as using the geodesic distance in eq. (13). One commonly used approach to determine the local distortion, is the Trustworthiness proposed in [31]. This measure provides a penalty when one of the $k$-nearest neighbors of an observation in the low-dimensional embedding is not one of the $k$-nearest neighbors in the original data space. For our particular case, we compare the low-dimensional embeddings with each other. Trustworthiness gives an index that goes from zero to one. The lower the distortion is, the closer Trustworthiness gets to one. Notice the trustworthiness does not need to be symmetric.

To measure global differences between embeddings, we employed the Mantel test [32], which gives a level of similarity between the distance matrices in the embedded dimensions. This gives an overall measure of the similarity between embeddings. The motivation for this test is the fact that distances are not independent from each other. For example, changing the position of a single observation will result in the distortion of $N-1$ distances. This implies that the direct calculation of the correlation between distances may not be enough to accurately assess the similarity between distance matrices. The Mantel test takes such dependencies between distances into account. The test outputs a correlation coefficient between 0 and 1, that can be interpreted similarly as the usual Pearson correlation.

4.3 Discussion

Figure 1AB shows how using distance (6), the embedding becomes robust with respect to $\gamma$ as $t$ increases. From a local perspective, the Trustworthiness measure among embeddings with different values of $\gamma$ gets closer to one using $t=10$ than for $t=1$. This shows a high local similarity independent of the value of gamma. From a global perspective, the Mantel test also shows a high dissimilarity for $t=1$, while for $t=10$ the embeddings become more similar. The previous assessments can be visually corroborated by looking at the embeddings themselves in Figure 1C. Finally, in Figure 1C we can observe a good representation of the sleep dynamics as there is a visually clear discrimination among sleep stages, especially for higher values of $\gamma$ or when using the diffusion geodesic distances. For these higher values of $\gamma$, the central structure of the embedding is more clearly defined than when using DM or lower values of $\gamma$. But largely, the embedding is robust to the choice of $\gamma$, suggesting that the specific choice of information distance is not too important when using the Mahalanobis distance in eq. (6) from the EIG framework. Furthermore, the results suggest that the use of the selected $f$-divergences does not provide any apparent advantage in this setting over the other information distances.

For the previous case, different values of $\gamma$ produced largely the same results. This is no longer the case when we use the Gaussian-based geodesic distance in eq. 8 as shown in Figure 2. First, Figure 2A and Figure 2B show quantitatively a greater difference between the embeddings than in Figure 1. This can be visually confirmed by looking at the embeddings in Figure 2C. In this case, the traditional DM tends to condense the structure together, and the use of the alternative $\gamma$ values may reveal more details of the structure of the data. The most left embedding is a clear representation of such a situation, where DM does not show a suitable discrimination of the sleep stages. But when increasing the value of gamma, a more suitable representation is achieved.

There are clear visual differences between the embeddings implementing distances (6) and (8). In Figure 3, we show a comparison of two embeddings colored by time steps. The left embedding is built using distance (6), achieving a good, clean visualization of the process across time. In fact, the different branches in the central region of the embedding are created from different periods of time, suggesting that transitions from different sleep stages may differ slightly depending on the total sleep time. In contrast, the Gaussian-based geodesic distance (8), displayed in the right embedding, not only inherits more of the original noise but also obscures the path of the process across time. Thus, using the Mahalanobis distance appears to better denoise the data and preserve the overall structure and time progression of this data.

5 Conclusion

In this work, we introduced a manifold learning tool for visualizing dynamical processes based on a diffusion framework. We addressed some of the shortcomings of the traditional diffusion maps approach for visualization, and used elements from PHATE and EIG to overcome them. We showed that when using the EIG-based distance, the visualization is robust to the choice of information distance.

We presented experimental results where we were able to discover sleep dynamics using solely EEG recordings. The 2-dimensional visualizations showed a clear distinction among sleep stages, as well as the time-varying progress of the processes.

Future work includes extending the analysis and comparison between different distance measures for time series data, further analyzing the impact of different information distances, and applying DIG to financial and biological data.