Extracting the Main Trend in a Data Set: The Sequencer Algorithm

As described above, we can achieve a meaningful ordering of the data set by minimizing the similarity or distance between adjacent objects. Doing so requires the choice of a distance measure. In order to be generic, we include several commonly used metrics. This will allow the algorithm to consider various aspects of the data and its features. Similarly, we consider a list of scales on which the relevant information can be distributed.

By default, we use the following metrics: (i) the Euclidean Distance, (ii) the Kullback-Leibler Divergence (KL Divergence; Kullback & Leibler 1951), (iii) the Monge-Wasserstein or Earth Mover Distance (EMD; Rubner et al. 1998), and (iv) the Energy Distance (Székely 2002). The definitions and properties of these metrics are described in Appendix A. If desired by the user, this list can be expanded to better suit a particular application. We note that this default set includes metrics with different properties: the Earth Mover Distance and Energy Distance are sensitive to the magnitude of displacements along the i coordinate. This is important with continuous measurements, for example, performed as a function of space or time, for which the derivative with respect to i carries relevant information. In contrast, the Euclidian Distance and the KL Divergence treat the different pixels of X_i as different dimensions and are insensitive to index shuffling. They provide a qualitatively different view of the information content.

The observable signature of an underlying trend can exist on different scales that may not be known a priori. In order to be generic, it is important to consider a range of scales. To do so, we perform a dyadic decomposition of each object X, i.e., we divide each one into a series of contiguous segments whose length is given by N_pix/2^l . This results in an ensemble of segments that allows us to look at each data object hierarchically, starting from its entirety (l = 0) and creating a binary tree such that the deepest scale corresponds to about 20 pixels. It is clear that one could sample the data using other strategies, for example, using specific convolution kernels. However, for simplicity and generality, the default version of the algorithm uses a dyadic decomposition. Thus, for a given metric k and scale l, the object is split into 2^l segments, and we refer to each segment using the index m. The maximum depth or the ways in which the data are decomposed can be modified by the user if need be.

Algorithm 1. Sequencer pseudo-code

set list of metrics

set list of scales

for each metric k do

for each scale l do

set list of segments

split objects X into m segments Xm

normalize each segment to have a sum of 1

for each segment m do

Dklm = distance matrix(Xm )

${\mathrm{MST}}_{\mathrm{klm}}$ = Minimum Spanning Tree(Dklm)

${\eta }_{\mathrm{klm}}={a}_{\mathrm{klm}}/{b}_{\mathrm{klm}}$ = elongation $({\mathrm{MST}}_{\mathrm{klm}})$

end

Dkl = η-weighted average of individual Dklm

${\mathrm{MST}}_{\mathrm{kl}}$ = Minimum Spanning Tree(Dkl)

${\eta }_{\mathrm{kl}}$ = elongation $({\mathrm{MST}}_{\mathrm{kl}})$

end

end

P = combined proximity matrix populated by η-weighted edges of MSTkl

for each pair of objects $i,j$ do

${D}_{{ij}}^{\mathrm{combined}}=1/{P}_{{ij}}^{\mathrm{combined}};$

end

$\mathrm{MST}$ = Minimum Spanning Tree(D)

Sequence = Breadth First Search path($\mathrm{MST}$)

Download table as:  ASCIITypeset images: 1 2

Our goal is to look at the data for each metric k, scale l, and segment m and estimate the level to which an underlying trend is present. To do so, we proceed as follows. First, for each scale l and each segment m, we first normalize each object such that the sum over its components is 1. By doing so, we only extract information from shape rather than amplitude. We then extract geometrical properties of the set of graphs characterizing similarities between all pairs of objects:

• 1.  
  Graph minimum spanning tree: for each metric k, scale l, and segment m, we compute an ${N}_{\mathrm{obj}}^{2}$ distance matrix D_klm, which represents a fully connected graph. We then compute its minimum spanning tree using Kruskal's algorithm (Kruskal 1956). It gives us a set of k × l × m trees with N_obj nodes. The key information on the presence of an underlying trend resides in the shape of these graphs.
• 2.  
  Graph length: The least connected node, j_LC, of the minimum spanning tree is expected to belong to its longest branch. To identify it, we compute the closeness centrality of each node (Freeman & Freeman 1978) and select the one with the smallest value. We then compute the shortest path Δ_klm (j_LC, j) between this node and every other node in the minimum spanning tree. The shortest path is a unitless integer counting the minimal number of edges between the two nodes (see the Appendix for additional details). We then define the major axis of the minimum spanning tree to be the average of the shortest paths over all nodes:

  $$\begin{eqnarray}&&{a}_{\mathrm{klm}}=\langle {{\rm{\Delta }}}_{\mathrm{klm}}({j}_{\mathrm{LC}},j){\rangle }_{\mathrm{node}\ j}.\end{eqnarray} \tag{ 1 }$$
• 3.  
  Graph width: Each node in the graph can be assigned to a level q that corresponds to a unique value of Δ_klm(j_LC, j), that is, the shortest path between all the nodes that are assigned to the level q and the least connected node: Δ_klm(j_LC, j) = Δ_q (see Appendix A.2 for an illustration). The width of a level q, ${{\rm{\Delta }}}_{\mathrm{klm}}^{\perp }(q)$, is defined as the number of nodes that are assigned to it. We use the average of this quantity as an estimate of the average half-width, or minor axis b, of the minimum spanning tree:

  $$\begin{eqnarray}&&{b}_{\mathrm{klm}}=\frac{1}{2}\, \langle {{\rm{\Delta }}}_{{klm}}^{\perp }(q){ \rangle }_{\mathrm{level}\ q}.\end{eqnarray} \tag{ 2 }$$
• 4.  
  Graph elongation: for each metric k, scale l, and segment m, we then define the elongation of the minimum spanning tree as its average height divided by its average width:

  $$\begin{eqnarray}&&{\eta }_{\mathrm{klm}}=\mathrm{elongation}({D}_{\mathrm{klm}})=\displaystyle \frac{{a}_{\mathrm{klm}}}{{b}_{\mathrm{klm}}}.\end{eqnarray} \tag{ 3 }$$

  This quantity can then be used to characterize the level to which the signature of a continuous trend is apparent in a given segment of the objects, through a given metric and on a given scale. Importantly, by being a ratio of numbers of edges, this parameter can be defined irrespective of the metric used. It is a summary statistics describing geometrical properties of the minimum spanning tree and, as a result, topological properties of the data.
• 5.  
  Aggregation of scales and metrics: each minimum spanning tree represents a sequence viewed through a given metric k and scale l of a segment of the data m. The elongation η_klm of its corresponding minimum spanning tree carries information on the level at which an underlying trend is detected. We can first combine the information obtained for all segments by creating a global distance matrix D_kl using an elongation-weighted average of our set of minimum spanning trees:

  $$\begin{eqnarray}&&{D}_{\mathrm{kl}}={\left\langle {\eta }_{\mathrm{klm}}.{D}_{\mathrm{klm}}\right\rangle }_{m}.\end{eqnarray} \tag{ 4 }$$

  This provides us with k × l different "views" of the data, which we can attempt to aggregate. Here we need to keep in mind that different metrics k will result in distance matrices D_kl with different units. To meaningfully combine information obtained from different metrics, we will only extract the topological information of the resulting minimum spanning trees, as given by their edge counts, and use an elongation-weighted average to create a "proximity" matrix:

  $$\begin{eqnarray}&&{P}^{\mathrm{combined}}={\left\langle {\eta }_{{kl}}.\#\,\mathrm{of}\,\mathrm{edges}(\mathrm{MST}({D}_{{kl}}))\right\rangle }_{{kl}}.\end{eqnarray} \tag{ 5 }$$

  We set all the elements that are not populated by the edges of the minimum spanning trees to zero (no connection between the corresponding nodes), and we set the elements on the diagonal to infinity (the proximity of a node to itself is infinite). This then allows us to define a combined distance matrix whose elements i, j are defined as ${D}_{{ij}}^{\mathrm{combined}}=1/{P}_{{ij}}^{\mathrm{combined}}$. This distance matrix D^combined provides us with a multiscale and multimetric characterization of the data set. We then compute its minimum spanning tree and corresponding elongation. Equation (5) performs an aggregation using a weighted average. We have found this approach to be appropriate in practice. One can of course consider examples for which most scales and metrics only contribute noise and dominate the final average. In such cases, the use of a max() would provide better results than a global average.
• 6.  
  Extraction of a final sequence: in order to extract a sequence from the minimum spanning tree of the combined distance matrix, we must select a particular walk within the tree, i.e., we must select the relative order in which we visit all the nodes within the graph. As done above, we define the starting point of the sequence using the least connected node of the combined minimum spanning tree. From this starting point, we walk through the graph using a Breadth First Search (BFS; Cormen et al. 2009). If the combined minimum spanning tree presents an appreciable elongation, a BFS traversal is expected to define the main trend in the data set. Thus, we expect the main branch of the tree to represent the sequence and secondary branches to represent the scatter or the secondary sequence.

We point out that the addition of one or several new objects to the data set can be done without redoing the full process. Once a sequence has been obtained, one can easily insert a new object into it, by performing a straightforward neighbor search. More details are provided in the Appendix.

The algorithm described above requires distance matrix calculations that scale as ${ \mathcal O }({N}_{\mathrm{obj}}^{2})$. For large data sets this approach can be computationally demanding. A useful alternative is to apply the technique to a subset of the data with N_s ≪ N objects, build the skeleton of a sequence, and then populate it with the remaining points. This process leads to a faster computation that allows one to process significantly larger data sets, at the cost of obtaining only an approximate result. A more detailed description of this method is provided in the Appendix.

We first construct a synthetic data set with a well-defined trend that exists only on small scales, on top of a varying background. We create 200 one-dimensional objects with N_pix = 400 presenting four narrow Gaussian pulses whose positions vary continuously from one object to another to form a clear sequence. To this we add random large-scale fluctuations using a Gaussian process. We show the shuffled data set in the top panels of Figure 2, where the left panel shows a subset of the objects in the sample and the right panel visualizes the full data set. Each row represents a different object, and the color-coding represents the relative intensity in each of its pixels. As can be seen, it is visually difficult to identify an underlying trend in this collection of objects. We apply the Sequencer to this data set and show its output in the bottom left panel. The algorithm is capable of detecting the overall structure in the data set even though only a small fraction of the pixels carry relevant information.

Zoom In Zoom Out Reset image size

For comparison, we show how t-SNE and UMAP, two of the most popular dimensionality reduction techniques, perform on the same simulated data set. To obtain a sequence using t-SNE or UMAP, we apply these techniques to embed the input data set into one dimension and then rank-order the objects according to their assigned value in this dimension. We point out that, as mentioned in the introduction, using these algorithms requires setting parameters. In practice, the user typically runs such an algorithm multiple times, varying those parameters until the "best" result is obtained. Here we present only the best sequences obtained after optimizing the distance metric in an automatic manner, using the elongation of the corresponding graphs. The details of this optimization are given in the Appendix. Clearly, t-SNE and UMAP fail to detect the sequence in narrow pulse locations, as they only consider an object as a whole, without the ability to focus on individual segments and understand that a well-defined trend exists on small scales.

To illustrate how the Sequencer identifies the scale and metric revealing a meaningful trend, in Figure 3 we present the resulting ordering of the algorithm for each metric and a set of scales as a function of the corresponding elongation parameter. The most elongated minimum spanning trees are obtained using the Euclidean Distance and KL Divergence measured over small scales (l = 4 and l = 5). These scales and metrics dominate the weighted averages in Equations (4) and (5) and, as a result, define the final ordering of the data.

Zoom In Zoom Out Reset image size

We now present a series of examples using one-dimensional data. For each data set, we first display 10 objects to convey the typical level of complexity. We then show a visualization of the entire randomly ordered data set, followed by the same data this time ordered by the Sequencer. We then display several objects from the ordered data set. In each case, we indicate the metrics and scales that were identified as the most informative, based on their contributions to the minimum spanning tree elongation-weighted averages in Equations (4) and (5).

Spectroscopic data. We start with a sample of 1000 spectra of stars from the publicly available Sloan Digital Sky Survey (York et al. 2000). Each spectrum is a measurement of the brightness of a star as a function of wavelength. The data set ordered by the Sequencer displays visible trends in both large-scale and small-scale features. A physical interpretation reveals that these continuous variations correspond to a sequence in the temperature of the stars. The top of the sequence is dominated by hot stars, which exhibit absorption lines due to hydrogen atoms, and the bottom part is dominated by cooler stars that exhibit absorption lines due to other, heavier elements.

We repeat the analysis with a set of 1000 spectra of quasars spanning a range of properties. This time we detect a trend in distance (redshift): the bright and narrow features shift continuously throughout the sequence. We point out that, using their default settings, most dimensionality reduction techniques (e.g., PCA, t-SNE, and UMAP) are insensitive to pixel shuffling and often fail to detect horizontal shifts in the data (see Appendix A.1).

Our third example shows another set of 1000 quasar spectra (Trump et al. 2006). The dark pixels correspond to flux deficits due to absorption by gaseous clouds present in front of the background light sources. We point out that this data set presents a higher level of apparent stochasticity than the previous example and is more difficult to interpret. Here, after reordering the data, we also smooth them in the Y direction, using a running median filter. This step compensates for the noise and makes weaker trends more easily apparent. We can observe dark regions corresponding to the absorption of light. The algorithm reveals the existence of two distinct populations of absorbers at the top and bottom of the ordered data set. This is most obvious when examining the data near λ = 1500 Å. These systems are known to exhibit different physical properties (e.g., Gibson et al. 2009 and references therein). This example illustrates how the algorithm can naturally perform a clustering task (and define sequences within each cluster), even in the presence of a substantial amount of noise.

Images from the natural world. In the case of simple, highly symmetric physical objects, it is sometimes possible to use a (physical) model and describe each object using a set of parameters. In such a case, it might be possible to identify an underlying trend in the data by looking at a trend in the best-fit parameters. When complexity increases, quantitative models based on a small number of parameters are no longer available. The search for trends requires a data-driven approach. To illustrate the performance of the algorithm in a higher-complexity regime where physical modeling of the data is out of reach, we use a set of images from the natural world (natural images). Due to the great variety of shapes and textures, the structural information of such images is expected to be distributed over a wide range of scales and cannot be described by a generic model. To stay in the simpler regime of one-dimensional objects, we randomly shuffle the rows of images and attempt to recover the original ordering. One example is shown in Figure 4. The left panel illustrates the complexity of each object and the lack of a simple model to characterize them. For the two examples shown, the algorithm is able to recover the original input. In the Appendix we present more examples using natural images and, for comparison, show the corresponding outputs for t-SNE and UMAP. We point out that, for representational purposes, all reordered images are flipped so that they match the expected orientation, when necessary. The algorithm has obviously no information regarding the correct orientation of the output.

Zoom In Zoom Out Reset image size

Displaying the spatial variation of vectors. Displaying spatial information, for example on a map, is often done through the use of a color bar, i.e., a sequence. This allows one to display a scalar value as a function of position. If the dynamic range is large, then techniques such as histogram equalization can be used to map the collection of values onto the appropriate range provided by the color bar. If the quantity to visualize is not a scalar but a vector or a distribution, there is no systematic way to define a mapping onto a color bar. However, by attempting to order these objects into a sequence, our algorithm provides us with a generic way to meaningfully assign a color to each vector or distribution. The sequencer does to vectors what histogram equalization does to scalars.

To illustrate how the algorithm can be used to display a set of distributions on a map, we apply it to a data set from structural seismology. Using the surface wave phase velocity at each of about 14,000 locations across the contiguous United States (Olugboji et al. 2017), we extract the Love wave dispersion curves specifying velocities at a set of 11 periods, which are primarily sensitive to crustal structure. To homogenize the data, we convert the velocity at each period to its standard score. The ordered-index obtained by the Sequencer can be color-coded. Using the same color bar, we can create a map that naturally shows the geographic patterns in crustal structure, shown in Figure 5.

Zoom In Zoom Out Reset image size

Prominent regions of thick sediment, such as the Williston and Denver basins and the Mississippi Embayment, are traced out at one end of the sequence. Regions with very contrasting crustal structure, such as the Sierra Nevada, High Rockies, and the Snake River Plain, are traced out at the other end of the sequence. We can also point out the correspondence between the geographic patterns revealed by the Sequencer and physiographic provinces identified at the surface (Fenneman 1928), shown by black lines in Figure 5. This example illustrates how the proposed algorithm can be used for mapping properties encoded as distributions rather than scalar values.

As mentioned above, if two types of objects are present in a data set, the algorithm can identify two distinct clusters and will present them sequentially in the final ordering of the data. This will also apply to cases for which a minority of objects, typically labeled as outliers, differ from the rest of the data. In other words, outliers are typically found at one end of the sequence, as, by definition, they present a large dissimilarity to the rest of the objects. One can also imagine some objects for which only a small fraction of the pixels differ from the rest of the population. Such objects or set of pixels would also be labeled as outliers. Interestingly, our approach allows us to detect them in a simple manner. Having an ordered sequence in hand, one can meaningfully perform an averaging/smoothing operation along that dimension and reduce the amount of fluctuations not related to the detected sequence. This has the advantage of enabling the detection of weak trends that are not easily visible in the randomly ordered data set. Then, once a smooth sequence has been defined, one can look at the pixel-level differences between the original sequence and the smooth counterpart. The statistics of these residuals can reveal objects for which a subset of the pixels differ from the expected underlying trend.

A number of considerations must be kept in mind when using the Sequencer: First, the ordering operation performed by the algorithm makes use of one definition of simplicity, based on the elongation of a minimum spanning tree (Equation (3)). While this always provides a well-defined ordering, it might not necessarily lead to the trend expected from model-based considerations, which are often based on a number of assumptions regarding features, scales, metrics, etc. If prior knowledge is available, it is advised to consider limiting the data to the region(s) where a meaningful variation is expected. Similarly, in some cases, rescaling the values of the input data might lead to better results, especially in cases for which the data present a high dynamic range.

The algorithm attempts to make use of information on different scales that, by default, are logarithmically spaced. While this approach is meant to be generic, it is possible to imagine cases for which this hierarchical decomposition might not be optimal. Segmenting the data using different strategies (e.g., using wavelets such that Fourier frequency windows do not overlap) might lead to better results. The user can modify the default sampling strategy if need be.

The algorithm attempts to find a trend using all fluctuations present in the data. These fluctuations can be due to useful structural information or due to random noise. Without a model, the algorithm cannot distinguish between them. Noise fluctuations are partially reduced when averaging is performed across segments of the data and scales (Equation (4)) but do present a limitation in identifying the underlying trend. If the data set presents a range of noise levels, it is advised to first apply the algorithm to a subset of the data with a higher signal-to-noise ratio. As described above, obtaining an ordered sequence offers another opportunity to perform an averaging operation along the sequence and further improve the characterization of the underlying trend.

In order to assess the level of similarity between pairs of objects, the algorithm uses a default set of four distance metrics. We briefly describe them here. This set can be expanded by the user if need be. As generically done, we will treat a pair of data set objects as two probability density distributions p(i) and q(i). An important point to consider is whether the index i encodes a dimension—and can therefore be shuffled, or the value of a coordinate, such as distance, time, energy, etc.—in which case the derivative with respect to i carries valuable information. Distance metrics insensitive to shuffling will not take this into account. When the index i encodes the value of a coordinate, for simplicity we assume that all objects in the set are sampled in the same manner. If not, interpolation/resampling techniques need to be applied prior to using the Sequencer algorithm. The four metrics used by default are as follows:

• 1.  
  The Euclidean Distance ( L₂ ), i.e., the familiar distance between two points in an Euclidean space, is the default metric used in many fields. In particular, it is the default metric used by the dimensionality reduction algorithms t-SNE and UMAP. Given two objects that are represented by the vectors p(i) and q(i), L₂ is given by

  $$\begin{eqnarray}&&{{\rm{L}}}_{2}(p,q)=\sqrt{\sum _{i}{\left[p(i)-q(i)\right]}^{2}}.\end{eqnarray} \tag{ 6 }$$

  L₂ is nonnegative, is symmetric, and equals zero only when the two vectors are identical. Furthermore, L₂ is insensitive to the relative order of the values in the vector and therefore does not use any (useful) information from the derivatives with respect to index i.
• 2.  
  The Kullback-Leibler Divergence (KL Divergence), also called the relative entropy, quantifies the degree of surprise involved in seeing one distribution, given another one. It is widely used in machine learning to compare the similarity between two objects. For two discrete random variables with PDFs p(i) and q(i), it is given by

  $$\begin{eqnarray}&&\mathrm{KL}(p| | q)=\displaystyle \sum _{i}p(i)\,\mathrm{log}\frac{p(i)}{q(i)}.\end{eqnarray} \tag{ 7 }$$

  The KL Divergence is nonnegative and equals zero only when the two functions are similar in every bin i. One can see that it is asymmetric, namely, KL(p∣∣q) ≠ KL(q∣∣p), and it reaches infinity if, in at least one bin i, q(i) = 0. We further note that the KL Divergence is not sensitive to the relative order of the bins.
• 3.  
  The Monge-Wasserstein or Earth Mover Distance (EMD) measures a distance between two distributions from an optimal transport point of view. Intuitively, if the distributions are interpreted as two different ways of piling up a given amount of dirt, the EMD is the minimum work required to turn one pile of dirt into the other, where work is the amount (mass) of dirt moved times the distance by which it moved. The EMD is widely used in computer vision, and specifically in content-based image retrieval, as it resembles remarkably well a human's visual perception (see, e.g., Rubner et al. 2000). Interestingly, for one-dimensional distributions, the solution of the optimal transport can be obtained simply by (Ramdas et al. 2015)

  $$\begin{eqnarray}&&\mathrm{EMD}(p,q)={\int }_{-\infty }^{\infty }| P(i)-Q(i)| {di},\end{eqnarray} \tag{ 8 }$$

  where P(i) and Q(i) are the cumulative distribution functions of the random variables p(i) and q(i). For higher-dimensional distributions, computing the EMD is more computationally demanding. The EMD is nonnegative and symmetric. In contrast to L₂ and the KL Divergence, the EMD is sensitive to the order and the value of the indices. It can therefore track and quantify translations and/or derivatives with respect to i. As a result, the EMD can capture certain aspects of similarities between objects to which L₂ and the KL Divergence might not be sensitive.
• 4.  
  The Energy Distance (ED) provides another measure of statistical distance between two probability distributions. Székely (2002) showed that for one-dimensional real-valued random variables, p(i) and q(i), with cumulative distribution functions P(i) and Q(i), it is equivalent to

  $$\begin{eqnarray}&&\mathrm{ED}(p,q)={\left(2{\int }_{-\infty }^{\infty }| P(i)-Q(i)| \right)}^{2}{di}.\end{eqnarray} \tag{ 9 }$$

  The ED is nonnegative and symmetric. Similarly to the EMD, it is sensitive to the order and values of the indices. The differences between the ED and the EMD are analogous to the differences between the L₁ and L₂ norms.

Here we provide a visual support to summarize the key quantities used by the algorithm and the corresponding terminology introduced in Section 3. Figure 6 shows the example of a minimum spanning tree, where the nodes are color-coded according to their centrality measure. The least connected node of this graph is indicated in the lower left corner of the figure. It provides the starting point to traverse the graph using a Breadth First Search walk, i.e., along the longest branch of the graph and scanning each branch along the way. The nodes are ordered in levels according to their distance (in units of edges) from the starting point (marked with numerical labels in the figure). These levels are used to estimate the average width and height of the tree, which are then used to estimate the elongation of the minimum spanning tree. The graph length (Equation (1)) and width (Equation (2)) are also illustrated in the figure.

Zoom In Zoom Out Reset image size

The algorithm described in Section 3 requires distance matrix calculations that scale as ${ \mathcal O }({N}_{}^{2})$. For large data sets this approach can become too computationally demanding. A useful alternative is to first apply the technique to a subset of the data, build the skeleton of a sequence, and then populate it with the remaining points. This process leads to a faster computation and enables the analysis of much larger data sets. It can be implemented as follows: for a data set with N objects, we first select a subset with N_s ≪ N objects for which the distance matrix calculations are computationally feasible. This provides us with a first sequence.

• 1.  
  Growing sequence: given a selected sequence, the algorithm selects a fraction f_A of objects distributed uniformly. We refer to them as anchor points.
• 2.  
  Adding new objects: to populate the growing sequence with a new object i, we first perform a low-resolution search where we measure the distance between this object and the N_A = f_A N_s anchor points and find the two nearest neighbors. We then perform a high-resolution search computing distances between object i and all the nodes from the evolving sequence located between the two nearest anchor points. These distances are populated into the proximity matrix, with a weight that corresponds to the aggregated minimum spanning tree elongation for a given scale (Equation (5)). This proximity matrix is converted into a distance matrix, after which a minimum spanning tree is constructed. We note that this part of the calculation involves very sparse matrices, as only a very small fraction of nodes and edges are considered. This process of populating the growing sequence can be done in parallel for a group of objects. Finally, the updated growing sequence is obtained by a Breadth First Search traversal of the minimum spanning tree. This concludes a single iteration of the population phase.
• 3.  
  Updated sequence: after having populated the growing sequence with a group of new objects, we obtain an updated sequence. We can then repeat the process with the remaining points. Here we point out that, at each iteration, the number of anchor points, defined to be a fraction f_A of the size of the coarse sequence, grows linearly with the size of the growing sequence. In order to ensure convergence, the number of objects that are populated in every iteration must be smaller than the number of anchor points.

This process avoids a full ${ \mathcal O }({N}_{}^{2})$ calculation and allows one to search for sequences in large data sets in a more efficient manner. This is done at the cost of obtaining only an approximate result for which the accuracy depends on the choice of initial subset size N_s and the fraction f_A of anchor points to use. With the example presented in Figure 5, it reduces the computing time by about two orders of magnitude. For more details regarding the implementation of this approximate and faster version of the code, see the companion Github repository: https://github.com/dalya/Sequencer/.

In order to compare the performance of the Sequencer to that of t-SNE or UMAP, we need to introduce a procedure to automatically select their parameters leading to the "best" sequence but without using the full machinery of the Sequencer. To do so, we proceed as follows: given a set of hyperparameters, each of these techniques can be applied to embed a collection of objects or vectors into a one-dimensional space and naturally obtained an ordered list. The question then is whether this list reveals a meaningful trend in the data set. To address this, it is unfortunately not possible to directly use the elongation parameter in that embedding, as, by definition, the corresponding manifold only has one dimension, so $\eta ^{\prime}$ is always 1. To meaningfully estimate a corresponding elongation parameter, we need to have access to more information. To do so, we use these embedding techniques to obtain projections in two dimensions. In the corresponding spaces, we can now meaningfully estimate the elongation of these manifolds and use this information as a figure of merit to assess the level to which a sequence is found. Here we point out that the elongation parameter estimated in this manner can serve as a figure of merit only when the two-dimensional embedding does not strongly depart from a one-dimensional manifold. It is meaningfully defined only when the two-dimensional distribution of objects can be simply represented by a major and minor axis. In cases involving clusters and/or outliers in the two-dimensional projection, the elongation parameter no longer informs on the quality of an expected sequence.

As described above, we can then use the elongation parameter as a figure of merit to find the metric and/or hyperparameters that produce the most elongated sequence or, in other words, the trend presenting the most elongated one-dimensional manifold, i.e., the highest degree of continuity. To illustrate this, we apply t-SNE and UMAP to two images for which the rows have been randomly shuffled. In Figure 7 we show how varying the hyperparameters of each technique affects the quality of the resulting embedding and how optimizing for the largest value of the elongation parameter naturally leads to the recovery of the original image. We also point out that optimizing the hyperparameters of these techniques does not always allow us to recover the original image.

Zoom In Zoom Out Reset image size

Having introduced the normalized elongation parameter to (i) automatically optimize the hyperparameters and/or metrics of t-SNE/UMAP and (ii) compare the resulting embeddings of different techniques, we now show that the Sequencer algorithm can identify an underlying trend in cases where both t-SNE and UMAP do not. This is similar to the example shown in Appendix B, but, this time, we are using real data rather than artificially generated distributions. To illustrate this point, we select images from the COCO data set (Lin et al. 2014). As done previously, we treat the rows as a collection of ordered objects whose order is randomly shuffled. We use them as inputs to the Sequencer, t-SNE, and UMAP, and in each case we attempt to recover the original ordering using the procedure described above. For t-SNE, we consider the four distance metrics and the hyperparameters: learning_rate = [10, 50, 100, 150] and perplexity = [10, 50, 70, 100]. For UMAP, we consider the same distance metrics and the hyperparameters: n_neighbor = [10, 50, 100, 150] and min_dist = [0.1, 0.2, 0.5, 0.8]. We note that this coarse sampling of the hyperparameters of these two techniques appears to be enough for this type of data. A higher-resolution optimization of the parameters does not provide better results.

The results are shown in Figure 8. For each input shuffled image, we present the output of the Sequencer, as well as the "best" outputs obtained with t-SNE and UMAP. In all cases, the elongation-based optimization leads to meaningful segments of the original images. Interestingly, in a number of cases, the Sequencer outperforms the best possible outputs obtained with t-SNE and UMAP. As described in Section 4.1, this is due to the ability of the Sequencer to look for a global trend using multiple metrics and a range of scales.

Zoom In Zoom Out Reset image size