---
abstract: 'When exploring large time-varying data sets, visual summaries are a useful tool to identify time intervals of interest for further consideration. A typical approach is to represent the data elements at each time step in a compact one-dimensional form or via a one-dimensional ordering. Such 1D representations can then be placed in temporal order along a time line. There are two main criteria to assess the quality of the resulting visual summary: *spatial quality* – how well does the 1D representation capture the structure of the data at each time step, and *stability* – how coherent are the 1D representations over consecutive time steps or temporal ranges? We focus on techniques that create such visual summaries using 1D orderings for entities moving in 2D. We introduce stable techniques based on well-established dimensionality-reduction techniques: Principle Component Analysis, Sammon mapping, and t-SNE. Our Stable Principal Component method is explicitly parametrized for stability, allowing a trade-off between the two quality criteria. We conduct computational experiments that compare our stable methods to various state-of-the-art approaches using a set of well-established quality metrics that capture the two main criteria. These experiments demonstrate that our stable algorithms outperform existing methods on stability, without sacrificing spatial quality or efficiency.'
author:
- |
    \
    [ ]{}
bibliography:
- 'references.bib'
---

Introduction
============

Time-varying data is ubiquitous and the size and the variety of the corresponding data sets is ever increasing. For the purpose of this paper, we distinguish three different types of time-varying data, which might appear in isolation or in combination. The first type are the quintessential *time-varying* data which reflect statistical attributes, such as population, or more general attributes such as size and price, in one or often many more dimensions. A second type are *dynamic* data which are characterized by insertions and deletions. The third type are *kinetic* data which have a spatial location (coordinates) that changes continuously over time and is sampled at discrete moments. Kinetic data typically represent moving entities. Any given time-varying data set can be both dynamic and kinetic.

![image](mrugsrelworklabeled.png){width="\linewidth"}

When exploring large time-varying data sets of any type, visual summaries are a useful tool to identify time intervals for further consideration. A typical approach is to represent the data elements at each time step in a compact form, which can then be placed in temporal order along a time line. For example, there are a variety of methods to handle time-varying graphs, such as Parallel Edge Splatting [@DBLP:journals/tvcg/BurchVBDW11] (Fig. \[fig:relworklabeled\]D) and Extended Massive Sequence Views [@van2014dynamic] (Fig. \[fig:relworklabeled\]E), that show the temporal evolution by drawing the graph at each time step in a narrow vertical strip. Similarly, Temporal Treemaps [@DBLP:journals/tvcg/KoppW19] (Fig. \[fig:relworklabeled\]F) encode hierarchies via (essentially) one-dimensional intervals and show the temporal evolution by placing these intervals consecutively along a line. Also Storyline Visualizations [@DBLP:journals/jgaa/DijkFFLMRSW17; @DBLP:journals/tvcg/LiuWWLL13] (Fig. \[fig:relworklabeled\]A) use a compact representation at each time step (essentially a pixel per protagonist); these representations must be coherent between consecutive time steps and as such trace a trajectory for each actor.

Arguably the most compact representation for one time step is a one-dimensional ordering of the data objects. Such an ordering directly translates to a grid-based visualization of associated attributes, where a vertical strip of $n$ grid cells encodes $n$ objects. Conceptually, computing a single linear order (linearization) for a given time step is a form of dimensionality reduction: after projecting to 1D, we sort the points by coordinate to obtain the order. Features visualized in the ordering enable an aggregated view on emerging trends, their composition and development over time. While the ordering prohibits the accurate perception of the spatial development of the tracked objects, a general spatial dynamic is encoded by the preservation of neighborhoods. Hence, spatial ordering approaches are best used for the in-situ visualization of large-scale kinetic data, enabling a quick estimation of general trends and developments. In principle, any aspect of the data can be used to create the ordering. For example, both MotionRugs [@BJC+19] (Fig. \[fig:teaser\]) and ParaGlide[@DBLP:journals/tvcg/BergnerSMAS13] (Fig. \[fig:relworklabeled\]B) compute orders from spatial locations for kinetic data, whereas Cui [*et al.*]{} [@DBLP:conf/apvis/CuiWLRMMG14] (Fig. \[fig:relworklabeled\]C) use node degree to order dynamic graph data.

For an ordering to meaningfully represent the input data, it needs to satisfy certain properties: (1) the ordering needs to accurately represent the data, as much as is reasonably possible; (2) the ordering needs to allow a user to follow data objects over time. For example, the ordering for dynamic graphs should capture the network structure, under insertions and deletions. For hierarchical data the ordering should capture the implied tree structure, under value changes and insertions and deletions. Finally, for kinetic data, the ordering should capture the spatial proximity under continuous movement of the entities. To facilitate the second requirement, the ordering needs to be temporally coherent and stable over consecutive time steps or temporal ranges. In this paper we focus on methods that create meaningful orderings for kinetic data.

[[**Formal problem statement.**]{}]{} Our input is a kinetic set $P = \{ p_1, \ldots, p_n\}$ of $n$ point objects moving in the plane. We sample their positions at $T$ consecutive time steps. That is, each object $p_i$ is a sequence of $T$ locations or points in the plane. We use $p_{i}(t)$ to denote the location of $p_i$ at time $t$, $1 \leq t \leq T$, and, correspondingly, $P(t)$ to denote the complete point set at time $t$. A visual summary $S$ of $P$ is a sequence of linear orders (permutations) of the points in $P$, one per time step. We denote the linear order at time $t$ by $S_t$. We consider $S_t$ to be a bijective function $S_t \colon P \rightarrow \{ 1, \ldots, n \}$. Thus, $S_t(p_i)$ denotes the rank of point object $p_i$ in the linear order at time $t$, and $S_t^{-1}(k)$ denotes the $k$^th^ point object in $S_t$. The quality of a visual summary $S$ is determined by two criteria:

Spatial quality.

:   How well does $S_t$ capture the spatial structure of $P(t)$? We characterize the spatial structure via local neighborhoods: we say that an order has high spatial quality if points that are spatially close in the input are also close in the order.

Stability.

:   How consistent are the linear orders over time? Here we can consider absolute changes between linear orders or changes in local neighborhoods, represented by nearest neighbors in the ordering. Both types of measures can be considered for consecutive time steps or over temporal ranges.

Clearly, a visual summary that uses the same order for all time steps is maximally stable. However, the spatial quality of this order will typically be low. Conversely, optimizing spatial quality for each time step in isolation tends to result in unstable summaries which make it more difficult for the user to track objects. Hence, we need algorithms that incorporate the temporal dimension into solving each time step. Ideally, techniques allow an explicit trade-off between spatial quality and stability of visual summaries.

[[**Contributions and organization.**]{}]{} We propose three dimensionality-reduction techniques to 1D that take stability into account. Two of our techniques require but simple modifications to existing methods: t-SNE and Sammon mapping. Our third novel algorithm, *Stable Principal Component* \[SPC\], is based on Principal Component Analysis (PCA) and allows an explicit trade-off between spatial quality and stability. The theoretical foundations for SPC lie in the stability analysis for kinetic shape descriptors by Meulemans, Verbeek and Wulms [@meulemans2019stability]. Despite the extensive underlying theory, SPC is simple and straightforward to implement. We describe all three methods in Section \[sec:stabledimreduc\] for kinetic data, but dynamic aspects (insertions and deletions) can be easily incorporated.

To order points, there are essentially three approaches via: (1) dimensionality reduction; (2) spatial subdivisions such as partitioning trees and space-filling curves; (3) hierarchical clustering. Generally, these methods are intended to perform well on a single time step and are not designed to be stable. In Section \[sec:algorithms\] we survey a representative set of state-of-the-art ordering methods.

To capture spatial quality and stability we rely on a set of well-established quality metrics. In particular, for spatial quality we use the so-called *Key Similarity* measures proposed by Guo and Gahegan [@guo2006spatial] that characterize spatial proximity via $k$ nearest neighbors. For stability we consider two different types of measures: absolute or neighborhood changes in the linear order. For absolute stability we use the number of *Jumps* and *Crossing* as proposed by Buchm[ü]{}ller [*et al.*]{} [@BJC+19]. We model neighborhood changes between orders via changes in the $k$ nearest neighbors, and again use the Key Similarity measures by Guo and Gahegan [@guo2006spatial]. We discuss metrics in detail in Section \[sec:metrics\].

In Section \[sec:experiments\] we report on computational experiments: we compare our stable algorithms with state-of-the-art ordering techniques, measuring the resulting quality using a variety of established metrics. Based on both real-world and synthetic data, these experiments demonstrates that our stable algorithms outperform existing methods on stability, without sacrificing spatial quality or computational efficiency. In Section \[sec:discussion\] we discuss the implications of our results, as well as current limitations and directions for future work.

Stable Dimensionality Reduction {#sec:stabledimreduc}
===============================

We describe three techniques for computing stable projections into a linear order for moving entities, which are based on well-established dimensionality-reduction techniques: PCA, Sammon mapping, and t-SNE. PCA has a clear indication of quality and representation, which permits interpolation between time steps. We use this fact to develop a stable PCA method which allows an explicit, user-configurable trade-off between spatial quality and stability (see Section \[sec:stablePCA\]). Sammon mapping and t-SNE are both based on gradient descent. In Section \[subsec:dimreduc\] we show how to leverage the existence of local minima to improve stability.

Stable Principal Component Analysis {#sec:stablePCA}
-----------------------------------

[r]{}[0.2]{} ![image](PCA.pdf){width="\linewidth"}

Our algorithm with an explicit, user-configurable trade-off between spatial quality and stability, uses the first principal component: an output of PCA originally introduced by Pearson [@pearson1901liii]. The first principal component is a vector in the direction along which the point set has most variance. Equivalently, it is the orientation of the line that minimizes the mean squared error to the points. As the most discriminative direction, we thus choose the first principal component vector to sort the points.

Meulemans, Verbeek and Wulms [@meulemans2019stability] study the stability of the first principal component from a theoretical point of view. Intuitively, they define stability as “small changes in the input lead to small changes in the output”. The paper analyzes the trade-off between spatial quality and the stability of various ways to summarize the shape of a point set. It shows that the first principal component exhibits unstable behavior, when the point set is not very stretched, and thus does not have a clear direction along which there is most variance. Our approach leverages this result by explicitly enforcing stability, whenever the variance along the first principal component is not high enough with respect to the variance in any other direction. The intuition behind this approach is as follows. If the variance along the first principal component vector is clearly higher than the variance along the second principal component, the direction is very discriminative: the point set is clearly stretched in this direction and sorting the points along this vector tends to lead to high spatial quality. If this is not the case, then the point set is “round” and the spatial quality is roughly equivalent for other directions as well. Our goal is to move between these directions of similar spatial quality in a stable way.

Meulemans, Verbeek and Wulms [@meulemans2019stability] introduce a “chasing algorithm”, which follows the optimal direction at fast as possible given a maximum rotation speed. This allows them to bound the solution quality and uses only the data of the “current time step”, which is particularly useful in a streaming setting. Here, however, we want to produce a visual summary and may reasonably assume that all data has already been recorded and is accessible. We can hence “look ahead” from the current time step to identify instabilities, before they actually occur.

[r]{}[0.2]{} ![image](SPC.pdf){width="\linewidth"}

[[**\[SPC$_\sigma$\] Stable Principal Component.**]{}]{} To create a stable version of principal component analysis, we use the optimal direction (first principal component) for any $t$ where $P(t)$ is stretched, as well as for the first and last time step. We then interpolate linearly between the orientations for these time steps to define directions also for any time steps that are not clearly stretched. Concretely, the Stable Principal Component algorithm is implemented as follows, using a parameter $\sigma$, $0 \leq \sigma \leq 1$, that determines when a point set is considered stretched. See Algorithm \[alg:stablepc\] for an overview.

The algorithm starts by finding the first principal component vector on $P(1)$ and then loops through time steps $t$, from $2$ to $T$ doing the following: We determine the first principal component vector $\textsc{pc}[t]$; as both $\textsc{pc}[t]$ and $-\textsc{pc}[t]$ describe the same orientation, we take $\textsc{pc}[t]$ to have the smallest absolute angle with respect to $\textsc{pc}[t-1]$, the direction of the previous time step. We add the signed angle between these vectors to $\alpha$, which accumulates the angle to interpolate over since the previous stretched (or first) time step. Let $v_1, v_2$ be the eigenvalues for the first and second principal component vectors; note that $v_2 \leq v_1$ by definition. Similar eigenvalues imply similar variances in those orthogonal directions: we determine how stretched $P(t)$ is via the ratio $r = \frac{v_2}{v_1}$. If $r > \sigma$, we consider the point set to be round and skip to the next time step. If $r \leq \sigma$ we consider it to be stretched. We remember the last time step $t'$ where the point set was stretched (or $t' = 1$). When we find a stretched point set at time $t$, we linearly interpolate the direction (first principal component) at time step $t'$ to the direction at time step $t$ for all skipped steps (round point sets), ensuring to turn the direction with an angle $\alpha$ in total (note that $\alpha$ can exceed $360^\circ$). After this loop, we have decided on a direction for each time step; we project the points at each time to the line defined by this direction and use that to infer the linear order $S_t$. The last time step $T$ is always considered “stretched” in the above, such that we end (as well as begin) with the actual first principal component.

Point set $P$ over $T$ time steps, and $\sigma \in [0,1]$ Visual summary $S$ for $P$

Set $\textsc{pc}[1]$ to the first principal component vector for $P(1)$ Set $t'$ to $1$ and $\alpha$ to $0$ Calculate first principal component vector $\textsc{pc}[t]$ and eigenvalues $v_1,v_2$ for $P(t)$ Add the signed angle between $\textsc{pc}[t]$ and $\textsc{pc}[t-1]$ to $\alpha$ Set $\textsc{pc}[t_s]$ to $\textsc{pc}[t']$ rotated over $\alpha \cdot \frac{t_s - t'}{t - t'}$ \[algline:interpolate\] Set $t'$ to $t$ and $\alpha$ to $0$ Define $S_t$ by projecting $P(t)$ on line through $\textsc{pc}[t]$ $S$

As Line \[algline:interpolate\] is executed at most once for each time step and computing two principal components and their eigenvalues of two-dimensional data takes linear time, the first loop runs in $O(nT)$ time. To sort the projections, the algorithm thus runs in $O(nT \log n)$ time, on $n$ points moving for $T$ time steps.

The explicit trade-off between spatial quality and stability can be configured via parameter $\sigma$. If $\sigma$ is set to a value close to $1$, the focus of the algorithm is on spatial quality, and only when the point set is very round, stability will be enforced; $\sigma=1$ eliminates interpolation and always uses the first principle component in every time step. However, if $\sigma$ is set closer to $0$, the focus will be on stability and even for moderately stretched point sets, linear interpolation can occur, thereby sacrificing spatial quality for stability; $\sigma = 0$ causes one interpolation, from the first principle component at $t = 0$ to the first principle component at $t = T$.

Gradient-descent methods {#subsec:dimreduc}
------------------------

The (spatial) quality of dimensionality reduction is a complex optimization function; as such, various local-search heuristics are commonly applied. For example, gradient descent is used by Sammon mapping to preserve distances, and by t-SNE to preserve local neighborhoods. We leverage the nature of gradient descent to increase stability and make it effective for creating visual summaries. There are other dimensionality-reduction techniques, such as MDS [@Kruskal1964] and Isomap [@Tenenbaum2000], but based on their optimization functions we believe that they give similar results. For example, in the Euclidean plane, classical MDS is equivalent to PCA. We first recall Sammon mapping and t-SNE for a single time step $t$, before explaining our adaptation for improved stability.

[r]{}[0.2]{} ![image](SAM.pdf){width="\linewidth"}

[[**\[SAM\] Sammon Mapping.**]{}]{} Sammon mapping [@Sammon1969] aims to preserve distances. Let $d_{ij}$ denote the Euclidean distance between points $p_i(t)$ and $p_j(t)$, denote the projected (1D) coordinates by $x_i$, and let $\delta_{ij} = |x_i - x_j|$. Sammon mapping computes coordinates $x_i$, attempting to minimize this cost function: $$C = \frac{1}{\sum_{1 \leq i < j \leq n} d_{ij}} \sum_{1 \leq i < j \leq n} \frac{(d_{ij} - \delta_{ij})^2}{d_{ij}}$$ We start with random initial coordinates and use a steepest (gradient) descent algorithm to minimize cost C.

[r]{}[0.2]{} ![image](SNE.pdf){width="\linewidth"}

[[**\[SNE\] t-Distributed Stochastic Neighbor Embedding.**]{}]{} The goal of t-SNE [@tsne2008] is to preserve local neighborhoods in the dimensionality reduction. Like Sammon mapping, it attempts to minimize a particular cost function. Again, let $d_{ij}$ denote the Euclidean distance between points $p_i(t)$ and $p_j(t)$. Similarities between points are captured by a probability distribution: $$\mathcal{P}_{j | i} = \frac{\exp\left(-\frac{d_{ij}^2}{2 \sigma_i^2}\right)}{\sum_{k \neq i} \exp\left(-\frac{d_{ik}^2}{2 \sigma_i^2}\right)}$$ The values $\sigma_i$ are chosen depending on the predefined *perplexity* $\kappa$ (see [@tsne2008] for details); in our experiments we use $\kappa = 40$. We further define $\mathcal{P}_{ij} = \frac{1}{2n}(\mathcal{P}_{j | i} + \mathcal{P}_{i | j})$ and we set $\mathcal{P}_{ij} = 0$ if $i = j$. Denote the projected (1D) coordinates by $x_i$, and define $\delta_{ij}$ as $$\delta_{ij} = \frac{(1 + |x_i - x_j|^2)^{-1}}{\sum_{k \neq l} (1 + |x_k - x_l|^2)^{-1}}$$ The cost function is defined by the Kullback-Leibler divergence as: $$C = \sum_{i \neq j} \mathcal{P}_{ij} \log\frac{\mathcal{P}_{ij}}{\delta_{ij}}$$ Finally, this cost function is minimized in the same way as for Sammon mapping: starting with random initial coordinates and using gradient descent[^1] to minimize cost C.

[[**Stability improvements.**]{}]{} Sammon mapping (SAM) and t-SNE (SNE) use random initial coordinates. To improve the stability of both algorithms we initialize them with the solution of the previous time step, resulting in two stable versions, **\[SAMp\]** and **\[SNEp\]**.

Recently, Rauber [*et al.*]{} [@DBLP:conf/vissym/RauberFT16] described Dynamic t-SNE: a more explicit way of making t-SNE stable over multiple time steps. Their approach performs a global optimization over all time steps simultaneously, using a separate copy of each point for each time step. They enforce temporal coherence by adding a term to the optimization function depending on the distance between two copies of the same point at consecutive time steps. For two reasons we were not able to include this algorithm in our experiments. First, it is very slow. The paper reports a running time of about 6 minutes per time step. Although a single time step of our data consists of only hundreds of points, we consider thousands of time steps, making their algorithm prohibitively slow for our experiments. Second, the implementation of Dynamic t-SNE rarely gives meaningful output when run on our data[^2], which is much more kinetic than the data experimented on in [@DBLP:conf/vissym/RauberFT16]. We further believe that Dynamic t-SNE would converge slowly on a time-varying data set with many time steps: it would take at least $T$ gradient descent iterations for frames that are $T$ time steps apart to affect each other. Since t-SNE is already known to converge quite slowly, the combination may simply require too many iterations to obtain a reasonable solution. Thus, Dynamic t-SNE exacerbates the usual downsides of t-SNE, namely black-box parameter tuning and slow convergence.

Other Techniques {#sec:algorithms}
================

In addition to dimensionality reduction, various other techniques can be found in literature that are used to compute linear orders for a set of points. We compare dimensionality-reduction techniques and our new adaptations not only to each other, but also to several of such existing algorithms. Most of these algorithms are not designed to be stable and typically consider different time steps in isolation. To select suitable algorithms for our comparison, we choose the algorithms that performed best in the experiments by Guo and Gahegan [@guo2006spatial] and Buchm[ü]{}ller [*et al.*]{} [@BJC+19]. We may classify these remaining algorithms based on how they compute a linear order: (1) via spatial subdivisions; (2) via clustering. Finally, we also include a baseline algorithm that is solely focused on stability. Figure \[fig:orderings\] shows an example of the orderings generated by a selection of the algorithms, including the dimensionality-reduction techniques, for one time step of our test data for reference. We give a short overview of the algorithm here, while a full explanation can be found in Appendix \[app:algorithms\].

![ Orderings using (SPC) and (SNE) dimensionality reduction, space-filling curves (HIL) and clustering (CLC).[]{data-label="fig:orderings"}](orderings2.png){width="1\linewidth"}

[[**\[FXD\] Fixed order.**]{}]{} This algorithm outputs the same arbitrary linear order for every time step and hence serves as reference baseline for our experiments. With FXD, each horizontal line always represents the same moving entity. [[**Spatial subdivisions.**]{}]{} Several well-known linearization approaches, which are primarily used for spatial-indexing applications, are based on iterating through some spatial subdivision. These approaches encompass tree data structures and space-filling curves. Many variations exist; see [@DBLP:journals/debu/LuO93] for an overview. Here we focus on four established, representative techniques from this area: \[HIL\] Hilbert curve [@hilbert1891ueber], \[ZOR\] Z-order curve, \[PQR\] Point Quadtree [@DBLP:journals/acta/FinkelB74; @DBLP:books/daglib/0032640] and \[RTR\] R-tree [@DBLP:conf/sigmod/Guttman84].

[[**Clustering.**]{}]{} Another approach is to first compute a hierarchical clustering on the point set, and then order the points in such a way that clusters stay together. These algorithms are defined by how the points are clustered, and how the linear order is computed from the clustering. We use \[CLC\] Complete Linkage Clustering [@gordon1987review] and \[SNN\] Shared Nearest Neighbors [@Jarvis1973] to cluster points and derive an ordering from the cluster hierarchy as follows.

The hierarchical clustering is represented by a tree with the individual points stored in the leaves. We aim to order to leaves of such a tree without changing the cluster structure, that is, by only changing the order of the children of any internal node. We follow the algorithm by Bar-Joseph [*et al.*]{} [@Bar-Joseph2003] to efficiently compute the optimal order that minimizes the length of the path formed by visiting the input points in that order.

Metrics {#sec:metrics}
=======

To evaluate the various methods for computing a visual summary, we need metrics to quantify spatial quality and stability.

Spatial Quality
---------------

Spatial quality measures the correspondence between $P(t)$ and the linear order $S_t$. We capture this by considering the local neighborhood of a point, as characterized by its nearest neighbors. One way to measure changes in local neighborhoods is using an evaluation of dimensionality reduction via persistent homology as introduced by Rieck and Leitte [@DBLP:journals/cgf/RieckL15]. However, we choose not to use this type of measure. While this approach is more recent than the measure we are using, it does not compare to older results, it is more complex, and most importantly it is an indirect approach. Hence, we use the *Keys Similarity* measures as described by Guo and Gahegan [@guo2006spatial] to directly measure the changes in nearest neighbors.

To simplify notation, we omit dependencies on time step $t$, as the metrics consider each time step in isolation. Thus, $P$ denotes a point set in the plane, and $S$ denotes a linear order. Let $n(i,j) \in P$ denote the $j^{th}$ nearest neighbor of $p_i$ in $P$, for each $j$ with $1 \leq j \leq k$ for some constant $k$. We use $r(i,j)$ to denote the neighbor rank in $S$ between $p_i$ and $n(i,j)$. However, the difference in rank $|S(n(i,j)) - S(p_i)|$ is not unique. There are two neighbors at rank difference 1, two at rank difference 2, until we reach one end of a linear order. To avoid arbitrariness, we do not break ties but rather consider each pair with the same rank difference to have the same value for $r(i,j)$. Thus, there are two nodes with $r(i,j) = 1$ (rank difference 1), two nodes with $r(i,j) = 3$ (rank difference 2), etc.

Generally, Keys Similarity at time $t$ is then defined as $$KS(P, S) = \frac{\sum_{p_i \in P} \sum^k_{j=1} w(i,j) \cdot r(i,j)}{\sum_{p_i \in P} \sum^k_{j=1} w(i,j)},$$ where $w(i,j)$ denotes the weight or importance of maintaining the $j^{th}$ nearest neighbor of $p_i$ at time $t$ – note that these weights need not be the same at every time step. We use two variants of Keys Similarity, as introduced by Guo and Gahegan [@guo2006spatial].

[[**\[KSra\] Rank-weighted Keys Similarity.**]{}]{} We define $w(i,j) = 1/j$ inversely proportional to the rank, such that maintaining the closest neighbors is considered more important than the more distant neighbors. This gives the following metric, where $H_k$ is the $k^{th}$ harmonic number: $$KSra(P,S) = \frac{\sum_{p_i \in P} \sum^k_{j=1} r(i,j) / j}{\sum_{p_i \in P} \sum^k_{j=1} 1/j} = \frac{\sum_{p_i \in P} \sum^k_{j=1} r(i,j) / j}{n \cdot H_k}$$

[[**\[KSdi\] Distance-weighted Keys Similarity.**]{}]{} We define $w(i,j) = 1/\| p_i - n_t(i,j) \|$ inversely proportional to the Euclidean distance, such that maintaining close neighbors is considered more important than distant neighbors. In contrast to KSra, this variant does not treat neighbors at (nearly) identical distances differently.

$$KSdi(P,S) = \frac{\sum_{p_i \in P} \sum^k_{j=1} r(i,j) / \| p_i - n_t(i,j) \|}{\sum_{p_i \in P} \sum^k_{j=1} 1/\| p_i - n_t(i,j) \|}$$

[[**Other facets.**]{}]{} Our metrics focus on combinatorial aspects of the position of the point objects. Spatial structure in general knows many other facets, such as distances and directions between points, as well as density. However, a linear order inherently does not lend itself to represent such concepts.

Stability
---------

Stability or temporal coherence measures the similarity between consecutive orders in $S$. In our evaluation, we use the following three measures for stability. The first two are based on absolute changes in the order and match the measures used by Buchmüller [*et al.*]{} [@BJC+19] to evaluate MotionRugs. The latter uses neighborhoods, based on the concepts by Guo and Gahegan [@guo2006spatial].

We aim to compare the similarity between two linear orders, $S_t$ and $S_{t+1}$ for each $t$ with $1 \leq t < T$. We could easily use the same metrics to compare nonconsecutive orders, but this provides little insight for such inherently sequential data. To consider the stability over a temporal range $[t,t']$, we use standard summary statistics (e.g., average, minimum, or maximum) over all consecutive pairs.

[[**\[JMP\] Jump distance.**]{}]{} We quantify the jump distance for a single point object $p_i$ as the difference between its ranks in the two orders, that is, $|S_t(p_i) - S_{t+1}(p_i)|$. The jump distance between two orders is then the sum over all jump distances for each point object. $$\mathit{JMP}_t(P, S) = \sum_{p_i \in P} |S_t(p_i) - S_{t+1}(p_i)|$$ The value for $\mathit{JMP}_t(P,S)$ lies between $0$ (perfectly stable) and $n(n-1)/2$ (complete inversion of the order).

[[**\[CRS\] Crossings.**]{}]{} Whereas JMP penalizes any change in the order, many points moving up together may not constitute much change. Instead we may count the number of inversions or crossings in the order, that is, the pairs $p_i$, $p_j$ for which $S_t(p_i) < S_t(p_j)$ and $S_t(p_i) > S_t(p_j)$. The metric $\mathit{CRS}_t(P,S)$ lies between $0$ (perfectly stable) and $n(n-1)/2$ (complete inversion of the order). Buchmüller et al. [@BJC+19] also use Kendall’s $\tau$ coefficient to evaluate stability. We choose to omit this, as it is equivalent to $1-2\cdot \mathit{CRS}_t(P,S)/(n(n-1)/2)$. That is, Kendall’s $\tau$ is the same as CRS up to normalization to the range $[-1,1]$.

[[**\[KSte\] Temporal Keys Similarity.**]{}]{} We may also take the same approach as for spatial similarity and consider the similarity of local neighborhoods in both orders. As distances are not inherently meaningful in the combinatorial order and simply correspond to ranking differences, we use only the rank-weighted version of Keys Similarity. Also for this metric $\mathit{KSte}_t(P,S)$, we do not break ties in either order, but rather give them the same rank.

Experimental Evaluation {#sec:experiments}
=======================

To demonstrate the effects of taking stability into account when creating visual summaries, we compare the three stable methods to the other ordering algorithms, presented in Sections \[sec:stabledimreduc\] and \[sec:algorithms\], in a quantitative evaluation. For each algorithm, we assess the spatial quality and stability of the computed visual summaries according to the metrics discussed in Section \[sec:metrics\]. The parameter for SPC is explored after settling on the most effective measures for spatial quality and stability. Furthermore, we report on run time measurements. Before discussing the results in detail, we describe our test data.

![Shoal movement over 80 seconds. Several behaviors can be observed: (straight) translation, turning, milling (circling).[]{data-label="fig:movement"}](overview.png){width="1\linewidth"}

Data {#sec:data}
----

For comparability, we use the same data as MotionRugs [@BJC+19] along with a synthetic data set generated using Netlogo [@netlogo]. The first data set tracks 151 fish of the Notemigonus crysoleucas species (Golden shiner). Golden shiner fish live in large groups called “shoals”, moving in coordination at almost any given time. The 151 fish were tracked optically while moving through a 2.1m by 1.2m shallow water tank, thus avoiding movement in the third dimension. The tank did not feature any obstacles or hindrances besides the side walls. Different movement patterns can be observed in the data, which allows us to test quality in different situations. Among these patterns are uniform group movements, partial and complete changes of direction, circular movement patterns and changes in group density, speed, and acceleration. 2000 frames of movement were recorded at a rate of 25 frames per second, resulting in 80 seconds of available collective movement data. For each frame, the spatial coordinates of each fish are recorded in a Cartesian coordinate system. During this period, the fish first move straight through the tank (Figure \[fig:movement\], 00:10), then turn, move downwards, turn again, move straight and finally enter a so-called milling formation, moving in a circular shape.

In addition to the full data set, we also zoom in on four excerpts (see Figure \[fig:strategycomparison\]). The first three excerpts emphasize typical movements of the entities: translation represents straight group movements with predominantly stable internal group structure, turns involve the group turning with some change in the group structure, and milling show the fish entering and maintaining circular movement patterns. The fourth excerpt does not feature a particular movement pattern but triggers so-called “phantom splits” [@BJC+19] for certain ordering methods, most notably HIL, PQR, or SNEp. The shoal of fish appears to split, but this is purely an artifact of the method and not reflecting the data.

The second data set is generated with Netlogo using the Flocking model [@netlogoflocking] from the openly available Models Library within the Netlogo application. Minimal adaptations were made to the model to ensure the boundaries of the canvas do not wrap around, and the trajectories of the moving entities could be extracted easily.

In this section we mainly focus on the fish data set, while highlighting results of the Netlogo data set only whenever these are distinctively different from the fish data set, which is during the parameter experiment. The full analysis for the Netlogo data set is given in Appendix \[app:netlogo\].

![Spatial-quality metrics: mean KSra (left) and KSdi (right) for all algorithms over all frames of the fish data set.[]{data-label="fig:spatialquality-chart"}](spatialquality-chart.pdf){width="\linewidth"}

![Stability metrics: mean JMP, CRS (left axis), and KSte (right axis) for all methods over all frames of the fish data set.[]{data-label="fig:stability-chart"}](stability-chart.pdf){width="\linewidth"}

![MotionRugs for all algorithms, color encodes speed (red = fast, blue = slow). Below each we show KSdi (yellow) and KSte (blue), capped at 37.5 and 6.33, respectively. Four highlighted sections indicate interesting excerpts. []{data-label="fig:strategycomparison"}](stratsoverview.png){width="0.99\linewidth"}

Running Time
------------

We implemented and executed all algorithms in Java 11 on a workstation with two Intel Xeon E5-2687W CPUs at 3.10GhZ, 16 Cores, 128GB Ram and an NVidia Quadro M600 GPU, running Windows 10. We measure the running time only for computing the orderings excluding reading input, color mapping and rendering. The running times range from a few milliseconds for the Z-Order curve (ZOR) to just over 8 hours for t-SNE (SNE). General observations include comparably good performance for the subdivision methods (ZOR, HIL, PQR, RTR), with values under one second. Only SPC variants are on par with this speed.

Quality Results
---------------

![SPC examples showing how $\sigma$ affects interpolation. Arrows roughly show which frames are interpolated. []{data-label="fig:spcoverview"}](spcoverview.png){width=".99\linewidth"}

Figures \[fig:strategycomparison\] and \[fig:spcoverview\] show MotionRugs for all algorithms for the complete fish data set. The MotionRugs are accompanied by a visualization of the mean KSdi and KSte values for each frame, cut off slightly above the mean values of most algorithms. This ensures that the differences between the average behavior of the algorithms becomes visible at a glance. Tables \[sup-tab:experiments\] and \[sup-tab:experiments2\] in Appendix \[app:summarystats\] provide summary statistics over all time steps and for each metric, for both data sets. Below, we first discuss spatial quality and stability statistics separately, along with a discussion of how the metrics are reflected in the full visual summaries and the four excerpts. We follow up with an exploration into the effects of the parameter value on the outcome of SPC and finally consider the trade-off between spatial quality and stability for all methods.

[[**Spatial quality.**]{}]{} Figure \[fig:spatialquality-chart\] compares the spatial-quality measures KSra and KSdi, as measured on all algorithms used in our experiments. For both measures lower values indicate higher spatial quality. Overall, we see that the KSra measurements are slightly lower for all algorithms, except SNEp where KSdi has a minimal edge over KSra. As expected FXD achieves the worst spatial quality. Furthermore, SNEp and the algorithms using spatial subdivisions are outperformed by the clustering algorithms, and other dimensionality-reduction techniques. Comparing the spatial quality of SPC to the algorithms that perform best on spatial quality, we see that SPC achieves comparable spatial quality. The choices for parameter $\sigma$ of SPC on the fish data set are 0, 1, and variables $a=0.35, b=0.53, c=0.78$. The choice for the intermediate values $a,b$ and $c$ is different for the two data sets and will be justified in the parameter exploration. Due to the strong correlation of both measures, we focus only on KSdi in the remainder.

[[**Stability.**]{}]{} Figure \[fig:stability-chart\] compares the stability measures: JMP, CRS and KSte. While JMP and CRS measure absolute changes between orders, KSte captures changes in local neighborhoods. For each measure lower values indicate higher stability. We see that CRS results in lower values than JMP, which is expected: two entities can jump to different positions in the next frame without crossing, but they cannot cross each other without jumping. We do see some differences between the two data sets, as opposed to the results for spatial quality. For FXD the result is again obvious: all measures are at their minimum. While JMP and CRS are generally low, CLC, SNN and SNE show very high numbers. Those three algorithms also perform worst according to the KSte metric. Another outlier that performs poorly on KSte is RTR, which also performs comparatively poorly on JMP and CRS. Of the remaining algorithms, the spatial subdivisions perform worst on KSte. The SAM, SAMp and SNEp algorithms and the SPC variants show similar and very low mean values of KSte. Again, we observe a strong correlation between the three metrics, and thus consider only KSte in the remainder.

[[**Excerpts.**]{}]{} Let us briefly consider the excerpts of motion (translation, split, turn, milling) as highlighted in Figure \[fig:strategycomparison\]. The ordering method clearly defines the resulting visual patterns. Some algorithms, such as CLC and SNE, produce visually rather cluttered results despite a good spatial quality. This clutter is a consequence of instability: the summaries fail to convey patterns over time despite individual frames being objectively good. For some algorithms, for example SAM or SPC with high parameter values, stability is challenged at key moments during the turn and milling patterns. The KSte values are very high for only a few frames. Here we clearly see that even when on average the stability is good, spikes in stability lead to visual artifacts that disrupt patterns that can arise over consecutive frames. Interestingly, while the spatial quality deteriorates for SNEp when it comes to the milling patterns, the other algorithms do not seem to be sensitive to that. Visually, the undesirable phantom split pattern [@BJC+19] can be identified in the subdivision methods (HIL, ZOR, PQR, RTR) and SNEp, while others do not seem to be prone to these kind of visual artifacts or generally produce too fuzzy visual results for such patterns to appear.

[[**SPC parameter.**]{}]{} We now investigate the parameter $\sigma$ of SPC and its effect on the results. We run SPC for values for $\sigma$, specifically, for $101$ different values from $0$ and $1$ with increments of $0.01$. As discussed before, we use KSdi to asses the spatial quality of the visual summaries, and specifically we use the mean over all frames. For stability we use the mean as well as the max KSte to quantify stability. As we saw before, mean KSte captures cohesion over time, while max KSte should be low to prevent visual artifacts from disrupting temporal patterns. The results for the fish data set are shown in Figures \[fig:spcoverview\] and \[fig:parameter\]. Note that the highest plotted value of $\sigma$ is $0.95$, while the lowest is $0.29$. Values above and below these extremes are identical to results with $0.95$ and $0.29$ respectively. The values of $\sigma$ that are indicated by labels in the figures are chosen as representatives and used in our other experiments.

![A comparison between the mean and mean (left) as well as max and mean (right) for KSte and for KSdi, for uniformly distributed parameter settings of SPC on the fish data set.[]{data-label="fig:parameter"}](parameter-mean.pdf){width="\linewidth"}

![A comparison between the mean and mean (left) as well as max and mean (right) for KSte and for KSdi, for uniformly distributed parameter settings of SPC on the fish data set.[]{data-label="fig:parameter"}](parameter-max.pdf){width="\linewidth"}

Overall, we see an inverse relation between stability and spatial quality. Values of $\sigma$ closer to $1$ result in better spatial quality, while values closer to $0$ sacrifice some spatial quality for more stability. This is to be expected, as SPC$_{1}$ always projects the fish to the first principal component; this will likely lead to the best spatial quality that can be achieved for any parameter value.

As $\sigma$ is decreased, SPC increasingly uses interpolated lines for projection instead. This interpolation smooths changes in angle of the line, but the projection reflects spatial relations less accurately as a result. When $\sigma$ drops below $0.30$, the interpolation happens purely between the first and last frame of the data set. Contrary to expectation, this negatively affects both spatial quality and stability: the first principal component rotates both clockwise and counterclockwise at varying speeds, not matching the uniform interpolation over such a long time period; as a result, the interpolated lines do not correspond at all to the first principal components, neither in angles nor in rotation direction. This mismatch in angles leads to poor spatial quality per frame, while the mismatch in rotation direction also decreases stability.

![A comparison between the max for KSte and the mean for KSdi for all algorithms on the fish data set.[]{data-label="fig:spatialqualityVSmaxstability-chart"}](spatialqualityVSstability-chart.pdf){width="\linewidth"}

![A comparison between the max for KSte and the mean for KSdi for all algorithms on the fish data set.[]{data-label="fig:spatialqualityVSmaxstability-chart"}](spatialqualityVSmaxstability-chart.pdf){width="\linewidth"}

The Netlogo data set shows similar results albeit more surprising. We again observe an inverse relation between stability and spatial quality. However, for $\sigma$ between $0.45$ and $0.59$ this relation is absent. Increasing $\sigma$ also leads to worse stability and worse spatial quality. During instabilities we observe that at certain frames where SPC projects to an interpolated line, the spatial quality is better than when we project to the first principal component. It is hence not unreasonable that interpolating less (and using first principal component more) when increasing $\sigma$ from $0.46$ to $0.56$ can negatively affect spatial quality. Below $0.40$ both spatial quality and stability change erratically: certain instabilities are no longer interpolated over, creating bigger intervals of consecutive interpolation.

[[**Trade-offs.**]{}]{} Our main goal is to investigate the trade-off between spatial quality and stability. Figure \[fig:spatialqualityVSstability-chart\] shows a scatterplot on the means of KDdi and KDte of all algorithms. Since lower values indicate better quality for both, methods in the bottom-left corner perform well on both aspects. In both figures SPC variants are colored in shades of red, Sammon mapping variants in blue and t-SNE variants in green. The most stable variants according to KDte have fully opaque colors, while the unstable variants have a lighter shade.

These results clearly show that methods based on spatial subdivisions (ZOR, HIL) and space-filling curves (PQR, RTR), albeit fast to compute, perform poorly on spatial quality and stability. The clustering methods (CLC and SNN) as well as SNE, on the other hand, perform quite well on spatial quality, but exhibit very poor stability. Recall that these methods are also slow to compute.

The fixed order (FXD) and SNEp are on the other extreme, having good stability, but very poor spatial quality. Furthermore, the strong influence of initialization for t-SNE stands out. When initialized with random coordinates (SNE), the spatial quality is very good, but the stability is extremely poor. On the other hand, initializing t-SNE with the embedding of the previous time step (SNEp) greatly improves stability, but spatial quality suffers greatly.

That leaves SAM, SAMp, and SPC variants, which perform well on both aspects. We note that SAM and SAMp perform very similarly on KDdi (difference of $0.03$), but SAMp performs significantly better in terms of stability. Finally, SPC variants also strike a good balance between spatial quality and stability. All SPC variants have slightly worse spatial quality than SAM variants, but improve stability. However, recall that SPC is significantly faster to compute than the Sammon mapping algorithms SAM and SAMp.

The overall composition remains similar, but differences in stability are highlighted. Note that SAM, SAMp and SNEp are deteriorating with respect to other methods; we can also see clear bursts of instability in Figure \[fig:strategycomparison\] for these methods.

Figure \[fig:spatialqualityVSmaxstability-chart\] also highlights stability differences between SPC variants. SPC$_1$ always uses the first principal component, which can behave erratically for round point sets, decreasing stability. Our new method SPC overcomes this problem by interpolating over these bursts of instability. Indeed, SPC is largely unaffected for lower parameter values, having the smallest standard deviation overall (see Table \[sup-tab:experiments\] in Appendix \[app:summarystats\]).

Overall, stable methods such as SAMp, SNEp, and SPC for parameter values lower than $1$, performs very well in terms of average and worst-case stability, while only marginally sacrificing spatial quality in the case of SAMp and SPC. SPC does so at a fraction of the computational cost necessary for more complex dimensionality reduction techniques. Considering all the above, we conclude that stable methods are the best for computing visual summaries for time-varying (kinetic) data.

Discussion & Future Work {#sec:discussion}
========================

The experiments indicate that our stable methods perform best, in particular SPC: it performs better in stability, and performs as well as or better than its competitors in terms of spatial quality and computational efficiency. One the one hand, we leverage interpolation in our adaptation of PCA to SPC, allowing explicit parametrization. SAM and SNE were modified by changing the initial state of the gradient-descent computation; while this generally improves stability, the effect on spatial quality depends on the sensitivity of the underlying measure to local minima.

[[**Movement characteristics.**]{}]{} Our data, by nature and design, describes moving entities that primarily form a single, mostly convex cluster. Thus, proximity is the primary concern for determining neighborhoods and hence indicates which entities should be close to each other in the linearization. By definition, the first principal component captures the most discriminating axis – for our single cluster data this is most indicative of neighborhoods, explaining the performance of our method in terms of spatial quality. Methods based on finding clusters suffer in quality (either spatially or temporally) as there are no clear clusters to exploit. With multiple clusters it may be desirable to separate the different clusters in the linear order. The first principal component may do that to some degree, but it is not difficult to imagine a situation where clusters project to the same interval along the principal component vector, interleaving those clusters in the order. By their nature, clustering-based methods will perform better in separating the clusters. But our experiments show that such methods will nonetheless struggle to find a good, stable order within the clusters.

Our method and cluster-based methods therefore seem complementary – the one suitable for ordering within a cluster and the other suitable for identifying higher-level structures. It may thus be of interest to develop a new hybrid algorithm by combining cluster detection with dimensionality reduction, as explored also in [@DBLP:journals/tvcg/WenskovitchCRHL18]. This algorithm should identify and maintain clear, well-separated clusters that also persist over time, and order points of individual clusters using our algorithm.

In the case of a complexly shaped cluster, we face yet another issue. Clustering detection might not be adequate to find the necessary structure. As above, a single, straight axis for determining the order does not necessarily capture proximity or neighborhood structure well and is hence likely to give unsatisfactory results as well. Perhaps methods from topological persistence can play an important role to identify the necessary structures in such clusters.

We leave the development and evaluation of such algorithms for more complex data as future work. Our results show the potential here, for adapting existing methods to explicitly consider stability.

[[**Beyond spatial data.**]{}]{} Our stable methods can be used in any situation with time-varying in at least two (numeric) dimensions, to determine the ordering. In a MotionRug, another dimension is then used to color the elements in each order. Such an approach may thus be useful for providing an overview also for abstract data. However, we expect it to be primarily useful when proximity (or more generally, neighborhoods) of items are meaningful in the dimensions used to derive principle components. Investigating precise conditions under which this approach is effective is left to future work.

[[**Overview-first.**]{}]{} Visual summaries are primarily an overview-first tool: an analyst is not expected to use only the summary. Visual summaries are intended to give an analyst a rough idea of what happens during the motion of the entities, as a first entry point to find time spans or sets of entities to further investigate. It is therefore important to understand how movement patterns relate to patterns visible in the summary and vice versa.

To ensure that *collective* movement of subgroups leads to observable patterns in a MotionRug, we need the attribute used for coloring to be similar for spatially close entities. In our data this is the case for speed and inherent in properties derived from the spatial arrangement. Without a relation between spatial proximity and attribute value, the colors in the MotionRug may jump and it becomes difficult to follow entities or subgroups. Furthermore, we may want to augment a MotionRug or a visual summary in general with information about its spatial and temporal quality. We have augmented our MotionRugs in Figure \[fig:strategycomparison\], using a simple bar chart to show spatial quality and stability per time step. Various other encodings could be considered, e.g. reducing the saturation of the colors or underlining the MotionRug with two lines where the pixel colors indicate the spatial and temporal quality. How to best visually convey the spatial and temporal quality, and how this effects user understanding are left to future work.

[[**Acknowledgments.**]{}]{} The authors wish to thank Prof. Dr. Iain Couzin and Dr. Alex Jordan of the Max-Planck-Institute for Ornithology, Radolfzell, Germany, and Dominik Jäckle for helpful insights and discussions, and for providing the fish data set.

Other Linearization Techniques {#app:algorithms}
==============================

In addition to dimensionality reduction, various other techniques can be found in literature that are used to compute linear orders for a set of points. We compare dimensionality-reduction techniques and our new adaptations not only to each other, but also to several of such existing algorithms. Most of these algorithms are not designed to be stable and typically consider different time steps in isolation. To select suitable algorithms for our comparison, we choose the algorithms that performed best in the experiments by Guo and Gahegan [@guo2006spatial] and Buchm[ü]{}ller [*et al.*]{} [@BJC+19]. We may classify these remaining algorithms based on how they compute a linear order: (1) via spatial subdivisions; (2) via clustering. Finally, we also include a baseline algorithm that is solely focused on stability. Figure \[app:orderings\] shows an example of the orderings generated by a selection of the algorithms, including the dimensionality-reduction techniques, for one time step of our test data for reference.

![ Orderings for one time step generated using linear (SPC) and nonlinear (SNE) dimensionality reduction, space-filling curves (HIL) and clustering (CLC).[]{data-label="app:orderings"}](orderings2.png){width="1\linewidth"}

[[**\[FXD\] Fixed Order.**]{}]{} This algorithm outputs the same arbitrary linear order for every time step and hence serves as reference baseline for our experiments. With FXD, each horizontal line always represents the same moving entity.

Spatial Subdivisions
--------------------

Several well-known linearization approaches, which are primarily used for spatial-indexing applications, are based on the principle of iterating through some spatial subdivision. These approaches encompass tree data structures and space-filling curves. We focus on four established, representative techniques from this area, though many variations exist; see [@DBLP:journals/debu/LuO93] for an overview.

[r]{}[0.2]{} ![image](HIL.pdf){width="\linewidth"}

![image](ZOR.pdf){width="\linewidth"}

[[**\[HIL\] Hilbert curve and \[ZOR\] Z-order curve.**]{}]{} The Hilbert curve [@hilbert1891ueber] is a continuous space-filling curve. It can be applied to cover a spatial region in arbitrary precision by repeating the construction pattern recursively. A set of points in space can then be linearized by sampling the curve and noting the order in which the points are encoded on the curve. Another representative of space-filling curves is the Z-order curve, which differs from the Hilbert curve in its geometrical construction pattern resembling a Z shape, where the space is partitioned in four quadrants in the order NW, NE, SW, SE (see figure on the right). Both approaches differ in neighborhood retention and construction complexity, as Lu and Ooi describe [@DBLP:journals/debu/LuO93]. In their comparison, Guo and Gahegan [@guo2006spatial] found that Hilbert curves avoid long jumps better than the Z-Order curve, which in turn outperforms the Hilbert curve in the average of the compared metrics. Since both produce visually different outcomes, we include both strategies in the comparison.

[r]{}[0.2]{} ![image](PQR.pdf){width="\linewidth"}

[[**\[PQR\] Point Quadtree.**]{}]{} Quadtrees [@DBLP:journals/acta/FinkelB74] partition space recursively in four parts, until each part contains only a single point. Consequently, sparse areas cause fewer splits than dense areas. Standard quadtrees divide the space in equal parts, while Point Quadtrees split at an input point and thus potentially unevenly in terms of area. To derive the 1D ordering, a depth-first tree-iteration strategy is used; given the neighborhood structure in the tree, this is more suitable than a breadth-first strategy. See [@DBLP:books/daglib/0032640] for details on tree-iteration strategies. The standard quadtree essentially reflects a Z-Order curve linearization if the same quadrant iteration is applied. Hence, we use the point quadtree variant which produces different orderings due to the intermittent partition.

[r]{}[0.2]{} ![image](RTR.pdf){width="\linewidth"}

[[**\[RTR\] R-tree.**]{}]{} In R-Trees [@DBLP:conf/sigmod/Guttman84] objects are stored recursively in minimum bounding rectangles (MBR). Each MBR can hold at most a predefined number of objects, thus ensuring a minimum fill. In comparison to quadtrees, more complex balancing is necessary, recomputing the MBRs, when the object limit is reached. Note that MBRs can overlap. Again, a depth-first iteration strategy is used to order points in an R-Tree.

Clustering
----------

Another method to compute a linear order from a point set is to first compute a hierarchical clustering on the point set, and then order the points in such a way that clusters stay together. Algorithms of this type are defined by two aspects: (1) how the points are clustered, and (2) how the linear order is computed from the clustering. In the algorithms we consider below, we always use the following method to compute the linear order from the clustering.

The hierarchical clustering is represented by a tree with the individual points stored in the leaves. We aim to order to leaves of such a tree without changing the cluster structure: we can change only the order of the children of any internal node. We follow the algorithm by Bar-Joseph [*et al.*]{} [@Bar-Joseph2003] to compute the order that minimizes the length of the path formed by visiting the input points in that order. The algorithm uses dynamic programming to efficiently find the optimal order for every subtree placing two specific points at the first and last position in the order.

[r]{}[0.2]{} ![image](CLC.pdf){width="\linewidth"}

[[**\[CLC\] Complete Linkage Clustering.**]{}]{} Initially, every point is considered as a separate cluster to be hierarchically merged in a bottom-up fashion [@gordon1987review]. We do so by repeatedly merging the closest two clusters, until we obtain a single cluster. Distance between clusters is measured as the distance between their farthest points.

![The two quality metrics: mean for KSra (left) and KSdi (right column) for all algorithms over all Netlogo data set frames.[]{data-label="fig:spatialquality-chart2"}](spatialquality-chart2.pdf){width="\linewidth"}

[r]{}[0.2]{} ![image](SNN.pdf){width="\linewidth"}

[[**\[SNN\] Shared Nearest Neighbors.**]{}]{} This clustering algorithm [@Jarvis1973] works the same as CLC, but it uses a different metric than Euclidean distance to measure the dis(similarity) between two points. For two points $p$ and $q$, we first count the number of points $x$ that are in the set of $k$ nearest neighbors for both $p$ and $q$. We then define the *shared nearest neighbor* (SNN) distance between $p$ and $q$ as $1/(x+1)$. The SNN clustering is computed using the SNN distance instead of the Euclidean distance. In case of ties in the SNN distance, we use the Euclidean distance to break ties. In our experiments we use $k = 10$.

Experimental evaluation Netlogo data set {#app:netlogo}
========================================

We go over the statistics for the Netlogo data set used in our experiments in more detail in this section. Again, we first consider the spatial quality and stability separately, followed by the parameter exploration for SPC. Finally we discuss the trade-off between spatial quality and stability as observed on the Netlogo data set.

[[**Spatial quality.**]{}]{} As can be seen in Figure \[fig:spatialquality-chart2\], while the absolute values for the Netlogo data set are higher than for the fish data set, the relative values are very similar. The fixed order gives very bad results on spatial quality, followed by SNEp and RTR. The spatial subdivision techniques all perform similarly, and are slightly better than the previously mentioned techniques. Of the remaining algorithms, the clustering techniques (CLC and SNN) perform slightly worse than all remaining dimensionality reduction techniques (SAM, SAMp, SNE and SPC). For the Netlogo data set we have chosen different parameter values for SPC, specifically $a=0.40, b=0.59$ and $c=0.62$.

![The three stability metrics: mean for JMP, CRS (left), and KSte (right) for all algorithms over all Netlogo data set frames.[]{data-label="fig:stability-chart2"}](stability-chart2.pdf){width="\linewidth"}

[[**Stability.**]{}]{} In Figure \[fig:stability-chart2\] we plot the stability statistics for the Netlogo data set. While the chart looks quite different from the stability chart for the fish data set, this is mostly due to the fact that JMP and CRS count the absolute number of changes in the orders, whereas KSte is normalized. Since the Netlogo data set behaves less stable than the fish data set, all metrics show higher values. However, the Netlogo data set also contains more moving entities, which increases the absolute number of changes even further.

Comparing the statistics of JMP and CRS, we see very similar performances of all algorithms, with SNE being the least stable, while SAMp is the most stable. On the fish data set SNE was the least stable method overall, and while it still has the highest number of absolute changes on the Netlogo data set, it performs a lot better according to KSte, meaning neighborhoods are perserved relatively well over time. Now we see that RTR performs worst on KSte, followed by the space-filling curves (HIL and ZOR) and the clustering algorithms (CLC and SNN) together with PQR. The SPC variants perform relatively worse than on the fish data set, with parameter values close to but lower than 1 being optimal for stability. Finally SAM, SAMp and SNEp are the most stable according to KSte.

![A comparison between the max for KSte and the mean for KSdi, for uniformly distributed $\sigma$ of SPC$_{\sigma}$ on the Netlogo data set.[]{data-label="fig:parameter-max2"}](parameter-mean2.pdf){width="\linewidth"}

![A comparison between the max for KSte and the mean for KSdi, for uniformly distributed $\sigma$ of SPC$_{\sigma}$ on the Netlogo data set.[]{data-label="fig:parameter-max2"}](parameter-max2.pdf){width="\linewidth"}

[[**Parameter experiment.**]{}]{} As already explained in the main text, the results of the parameter experiment are slightly different for the Netlogo data set. In Figures \[fig:parameter-mean2\] and \[fig:parameter-max2\] the results are plotted. For the Netlogo data set, the cut-off values are $0.76$ and $0.32$, so everything above $0.76$ uses exactly the first principal component per frame, and similarly for values below $0.32$ we always interpolate between the first and last frame. The parameter values that are indicated by labels in the figures are the values we used in our other experiments. Note that there are two blue labels, which represent other values of interest that we will use in the analysis below.

We will consider the results between the values indicated by black labels in the figures. Starting from the lowest parameter value $0.32$ we see that increasing $\sigma$ has chaotic effects on both spatial quality and stability up to $0.40$ where this fickle behaviour ends. On closer inspection, the values between $0.32$ and $0.40$ constantly pick up more frames where the entities are stretched enough to use the actual first principal component. Since the intervals between which interpolation happens, constantly change, the results do not steadily change, but are quite erratic. Parameter value $0.38$ shows the worst combination, having both bad spatial quality and stability (max and mean).

![A comparison between the max for KSte and the mean for KSdi for all algorithms on the Netlogo data set.[]{data-label="fig:spatialqualityVSmaxstability-chart2"}](spatialqualityVSstability-chart2.pdf){width="\linewidth"}

![A comparison between the max for KSte and the mean for KSdi for all algorithms on the Netlogo data set.[]{data-label="fig:spatialqualityVSmaxstability-chart2"}](spatialqualityVSmaxstability-chart2.pdf){width="\linewidth"}

From $0.40$ to $0.45$ we see a steady decrease in stability and increase in spatial quality, as expected when increasing $\sigma$. Further increasing the parameter to $0.59$ has negative effects on both the spatial quality and stability. On closer inspection, this increase in spatial quality can be attributed to properties of the Netlogo data set. During instabilities we can observe that at certain frames where SPC projects to an interpolated line, the spatial quality is better than when we project to the first principal component. While we expect projections to the first principal component to have high spatial quality, it is not always the case, as we see here. In the Netlogo data set this occurs when the cluster of points changes direction and shortly does not form a convex shape. It is therefore not unreasonable that interpolating less (and using first principal component more) when increasing $\sigma$ from $0.46$ to $0.56$ can negatively effect spatial quality.

At $0.61$ SPC splits the interpolation over the two consecutive instabilities that were seen as one big instability. This split improves both spatial quality and stability, up to $0.62$ where there are a couple of non-interpolation frames between the instabilities. Increasing the parameter further leads to certain instabilities not being interpolated over any longer, which negatively affects the maximum KSte values observed for those runs of SPC.

Summary Statistics {#app:summarystats}
==================

The two tables on the next page provide the statistics over all time steps, for each metric and on every algorithm in our experiment. Each table shows the results for one of the data sets, the fish data set and Netlogo data set respectively. The data in the tables is used throughout the paper in a variety of charts and diagrams.

In the main text we explained that t-SNE was implemented from scratch, using the simplest form of the algorithm. There are extensions that approximate the gradient during gradient descent to improve run times. These extensions are integrated in most libraries that are currently available. When using the libraries we indeed saw faster run times for both SNE and SNEp, but stability was influenced by the approximations used in the libraries. We therefore chose to implement t-SNE from scratch to show its true capability to produce stable visual summaries. While this lead to slower run times than one would expect from state-of-the-art t-SNE, the run times of the libraries still greatly exceeded the run times observed for the other algorithms in our experiments.

[^1]: We tried using the existing implementation at <https://github.com/lejon/T-SNE-Java> to compute the t-SNE mapping. This implementation uses approximations to speed up the computation, which lead to artifacts in our results. We therefore implemented the default version of t-SNE ourselves. See Appendix \[app:summarystats\] for more details.

[^2]: The implementation often gives NaN as output. The authors [@Telea2019] have verified that this is a known problem with the implementation.
