A Spectral Method for Assessing and Combining Multiple Data Visualizations

In order that the proposed method is invariant to the respective scale and coordinate basis (i.e., directionality) of the low-dimensional embeddings generated from different visualization method, we start by considering the normalized pairwise-distance matrix for each visualization.

Specifically, for each $k\in\{1,2,...,K\}$, we define the normalized pairwise-distance matrix

where

is the un-normalized Euclidean distance matrix, and $\mathbb{D}^{(k)}=\text{diag}(\|\mathbb{P}^{(k)}_{1.}\|_{2},\|\mathbb{P}^{(k)}_{2.}\|_{2},...,\|\mathbb{P}^{(k)}_{n.}\|_{2})$ is a diagonal matrix with its diagonal entries being the $\ell_{2}$-norms of the rows $\{\mathbb{P}^{(k)}_{1.},...,\mathbb{P}^{(k)}_{n.}\}$ of $\mathbb{P}^{(k)}$. As a result, the normalized distance matrix $\bar{\mathbb{P}}^{(k)}$ has its rows being unit vectors, and is invariant to any scaling and rotation of the visualization $\{\mathbb{X}_{i}^{(k)}\}_{1\leq i\leq n}$.

The normalized distance matrices $\{\bar{\mathbb{P}}^{(k)}\}_{1\leq k\leq K}$ summarize the candidate visualizations in a compact and efficient way. Their scale- and rotation-invariance properties are particularly useful for comparing visualizations produced by distinct methods.

Our spectral method for assessing multiple visualizations is based on the normalized distance matrices $\{\bar{\mathbb{P}}^{(k)}\}_{1\leq k\leq K}$. For each $i\in\{1,2,...,n\}$, we define the similarity matrix

which summarizes the pairwise similarity between the candidate visualizations with respect to sample $i$. By construction, the entries of $\mathbb{G}_{i}$ are inner-products between unit vectors, each representing the normalized distances associated with sample $i$ in a candidate visualization. Naturally, a larger entry $(\bar{\mathbb{P}}_{i.}^{(k_{1})})^{\top}\bar{\mathbb{P}}_{i.}^{(k_{2})}$ indicates higher concordance between the two candidate visualizations. Then, for each $i\in\{1,2,...,n\}$, we define the vector of eigenscores $\widehat{\mathbb{s}}_{i}=(\hat{s}_{i,1},...,\hat{s}_{i,K})$ for the candidate visualizations with respect to sample $i$ as the absolute value of the eigenvector $\widehat{\mathbb{u}}_{i}\in\mathbb{R}^{K}$ of $\mathbb{G}_{i}$ associated to its largest eigenvalue, that is,

where the absolute value function $|\cdot|$ is applied entrywise. As will be explained later (Section 2.3.2), the nonnegative components of $\widehat{\mathbb{s}}_{i}$ quantify the relative performance of $K$ candidate visualizations with respect to sample $i$, with higher eigenscores indicating better performance. Consequently, for each candidate visualization $\{\mathbb{X}_{i}^{(k)}\}_{1\leq i\leq n}$, one obtains a set of eigenscores $\{\hat{s}_{i,k}\}_{1\leq i\leq n}$ summarizing its performance relative to other candidate visualizations in a sample-wise manner. Ranking and selection among candidate visualizations can be achieved based on various summary statistics of the eigenscores, such as mean, median, or coefficient of variation, depending on the specific applications. In particular, when some candidate visualizations are produced by the same method but under different tuning parameters, the eigenscores can be used to select the most suitable tuning parameters for visualizing the dataset. However, a more substantial application of the eigenscores is to combine multiple data visualizations into a meta-visualization, which has improved signal-to-noise ratio and higher resolution of the structural information contained in the data.

Importantly, as will be shown later (Section 2.3.5), the eigenscores essentially take the underlying true signals rather than the noisy observations $\{\mathbb{Y}_{i}\}_{1\leq i\leq n}$ as its referential target for performance assessment, making the method easier to implement and less susceptible to the effect of noise in the original data (Section 2.3.5 and Supplement Figure 21).

Using the above eigenscores, one can construct a meta-distance matrix properly combining the information contained in each candidate visualization. Specifically, for each $i\in\{1,2,...,n\}$, we define the vector of meta-distances with respect to sample $i$ as the eigenscore-weighted average of all the normalized distances respect to sample $i$, that is,

Then, the meta-distance matrix is defined as $\bar{\mathbb{P}}^{m}\in\mathbb{R}^{n\times n}$ whose $i$-th row is $\bar{\mathbb{P}}_{i.}^{m}$. To obtain a meta-visualization, we take the meta-distance matrix $\bar{\mathbb{P}}^{m}$ and apply an existing visualization method that allows for the meta-distance $\bar{\mathbb{P}}^{m}$ (or its symmetrized version $\bar{\mathbb{P}}^{m}+(\bar{\mathbb{P}}^{m})^{\top}$) as its input.

Intuitively, for each $i=1,2,...,n$, we essentially apply a principal component (PC) analysis to the normalized distance matrix $\bar{\mathbb{P}}_{i}=\begin{bmatrix}\bar{\mathbb{P}}_{i.}^{(1)}&\bar{\mathbb{P}}_{i.}^{(2)}&...&\bar{\mathbb{P}}_{i.}^{(K)}\end{bmatrix}\in\mathbb{R}^{n\times K}.$ Specifically, by definition (Jolliffe and Cadima, 2016) the leading eigenvector $\widehat{\mathbb{u}}_{i}$ of $\mathbb{G}_{i}=\bar{\mathbb{P}}_{i}^{\top}\bar{\mathbb{P}}_{i}\in\mathbb{R}^{K\times K}$ is the first PC loadings of $\bar{\mathbb{P}}_{i}$, whereas the first PC is defined as the linear combination $\bar{\mathbb{P}}_{i}\widehat{\mathbb{u}}_{i}=\sum_{k=1}^{K}\hat{u}_{i,k}\bar{\mathbb{P}}_{i.}^{(k)}$. Under the condition that the first PC loadings are all nonnegative (which is ensured with high probability under condition (C2) below), the first PC $\bar{\mathbb{P}}_{i}\widehat{\mathbb{u}}_{i}$ is exactly the meta-distance $\bar{\mathbb{P}}_{i.}^{m}$ defined in (2.9) above. When interpreted as PC loadings, the leading eigenvector $\widehat{\mathbb{u}}_{i}$ of $\mathbb{G}_{i}$ contains weights for different vectors $\{\bar{\mathbb{P}}_{i.}^{(k)}\}_{1\leq k\leq K}$ so that the final linear combination $\bar{\mathbb{P}}_{i}\widehat{\mathbb{u}}_{i}$ has the largest variance, that is, summarizes the most information contained in $\bar{\mathbb{P}}_{i}$. It is in this sense that the meta-distance $\bar{\mathbb{P}}_{i.}^{m}$ is a consensus across $\{\bar{\mathbb{P}}_{i.}^{(k)}\}_{1\leq k\leq K}$.

For our own numerical studies (Section 2.2), we used UMAP for meta-visualizing datasets with cluster structures, and used kPCA for meta-visualizing all the other datasets with smoother manifold structures, such as trajectory, cycle, or mixed structures. The choice of UMAP in the former case was due to its advantage in treating large numbers of clusters without requiring prior knowledge about the number of clusters (McInnes et al., 2018; Cai and Ma, 2021); whereas the choice of kPCA in the latter case was rooted in its advantage in capturing nonlinear smooth manifold structures (Ding and Ma, 2022). In each case, the hyper-parameters used for generating the meta-visualization were determined without further tuning – for example, when using UMAP for meta-visualization, we set the hyper-parameters the same as those associated to the UMAP visualization which achieved higher median eigenscore than other UMAP visualizations. Moreover, while in general UMAP/kPCA works well as a default method for meta-visualization, our proposed algorithm is robust with respect to the choice of this final visualization method. In our numerical analysis (Section 2.2), we observed empirically that other methods such as t-SNE and PHATE could also lead to meta-visualizations with comparably substantial improvement over individual candidate visualizations in terms of the concordance with the underlying true low-dimensional structure of the data (see Supplement Figure 11). In addition, the meta-visualization shows robustness to potential outliers in the data (Figures 4 and 5).

In Section 2.3, under a generic signal-plus-noise model, we obtain explicit theoretical conditions under which the performance of the proposed spectral method is guaranteed. Specifically, we show the convergence of the eigenscores to a desirable concordance measure between the candidate visualizations and the underlying true pattern, characterized by their associated pairwise distance matrices of the samples. In addition, we show improved performance of the meta-visualization over the candidate visualizations in terms of their concordance with the underlying true pattern, and its robustness against possible adversarial candidate visualizations. These conditions provide proper interpretations and guidance on the application of the method, such as how to more effectively prepare the candidate visualizations (Section 2.3.6). For clarity, we summarize these technical conditions informally as follows, and relegate their precise statements to Section 2.3.

• (C1’)

  The performance of candidate visualizations are sufficiently diverse in terms of their individual distortions from the underlying true structures.

• (C2’)

  The candidate visualizations altogether contains sufficient amount of information about the underlying true structures.

Intuitively, Condition (C1’) concerns diversity of methods in producing candidate visualizations, whereas Condition (C2’) is related to the quality of the candidate visualizations. In practice, Condition (C2’) is satisfied when the signal-to-noise ratio in the data, as described by (2.11), is sufficiently large, so that the adopted visualization methods perform reasonably well on average. On the other hand, from Section 2.3, a sufficient condition for (C1’) is that, at most $\sqrt{K}$ out of $K$ candidate visualizations are very similarly distorted from the true patterns in terms of the normalized distances $\bar{\mathbb{P}}^{(k)}$. This would allow, for example, groups of up to 3 to 4 candidate visualizations out of 10 to 15 visualizations being produced by very similar procedures such as the same method under different hyper-parameters.

To demonstrate the wide range of applicability and the empirical advantage of the proposed method, we consider visualization of three families of noisy datasets, each containing a distinct low-dimensional structure as its underlying true signal. We assess performance of the eigenscores and the quality of the resulting meta-distance matrix based on 16 candidate visualizations produced by multiple visualization methods.

For a given sample size $n$, we generate $p$-dimensional noisy observations $\{\mathbb{Y}_{i}\}_{1\leq i\leq n}$ from the signal-plus-noise model $\mathbb{Y}_{i}=\mathbb{Y}_{i}^{*}+\mathbb{Z}_{i},$ where $\{\mathbb{Y}_{i}^{*}\}_{1\leq i\leq n}$ are the underlying noiseless samples (signals), and $\{\mathbb{Z}_{i}\}_{1\leq i\leq n}$ are the random noises. Specifically, we generate true signals $\{\mathbb{Y}_{i}^{*}\}_{1\leq i\leq n}$ from various low-dimensional structures isometrically embedded in the $p$-dimensional Euclidean space. Each of the low-dimensional structures lie in some $r$-dimensional linear subspace, and is subject to an arbitrary rotation in $\mathbb{R}^{p}$, so that these signals are generally $p$-dimensional vectors with dense (nonzero) coordinates. Then we generate ${i.i.d.}$ noise vector $\mathbb{Z}_{i}$ from the standard multivariate normal distribution $\mathcal{N}({\bf 0},{\bf I}_{p})$, and use the $p$-dimensional noisy vector $\mathbb{Y}_{i}=\mathbb{Y}_{i}^{*}+\mathbb{Z}_{i}$ as the final observed data. In this way, we simulated noisy observations $\{\mathbb{Y}_{i}\}_{1\leq i\leq n}$ of an intrinsically $r$-dimensional structure. For our simulations, for some given signal-to-noise ratio (SNR) parameter $\theta>0$, we generate $\{\mathbb{Y}^{*}_{i}\}_{1\leq i\leq n}$ uniformly from each of the following three structures:

• (i)

  Finite point mixture with $r=5$: $\{\mathbb{Y}^{*}_{i}\}_{1\leq i\leq n}$ are independently sampled from the discrete set $\{\bm{\gamma}_{1},\bm{\gamma}_{2},...,\bm{\gamma}_{r+1}\}\subset\mathbb{R}^{p}$ with equal probability, where $\bm{\gamma}_{i}$’s are arbitrary orthogonal vectors in $\mathbb{R}^{p}$ with the same length, i.e., $\|\bm{\gamma}_{i}\|_{2}=\theta$ for $1\leq i\leq r+1$.

• (ii)

  “Smiley face” with $r=2$: $\{\mathbb{Y}^{*}_{i}\}_{1\leq i\leq n}$ are generated independently and uniformly from a two-dimensional “smiley face” structure (Supplement Figure 6 left) of diameter $\theta$, isometrically embedded in $\mathbb{R}^{p}$ and subject to an arbitrary rotation.

• (iii)

  “Mammoth” manifold with $r=3$: $\{\mathbb{Y}^{*}_{i}\}_{1\leq i\leq n}$ are generated independently uniformly from a three-dimensional “mammoth” manifold (Supplement Figure 6 right) of diameter $\theta$, isometrically embedded in $\mathbb{R}^{p}$ and subject to an arbitrary rotation.

The thus generated datasets cover diverse structures including Gaussian mixture clusters (i), mixed-type nonlinear clusters (ii), and a connected smooth manifold (iii). As a result, the first family of datasets was set to have $p=500$ and $n=900$, and were obtained by fixing various values of the SNR parameter $\theta$, and generating $\mathbb{Y}^{*}_{i}\in\mathbb{R}^{p}$ from the above setting (i) to obtain the noisy dataset $\{\mathbb{Y}_{i}\}_{1\leq i\leq n}$ as described above. Similarly, the second and the third families of datasets were obtained by drawing $\mathbb{Y}^{*}_{i}\in\mathbb{R}^{p}$ from the above settings (ii) and (iii), respectively, and generating datasets $\{\mathbb{Y}_{i}\}_{1\leq i\leq n}$ with $p=300$ and $n=500$, for various values of $\theta$.

For each dataset $\{\mathbb{Y}_{i}\}_{1\leq i\leq n}$, we consider $12$ existing data visualization tools including principal component analysis (PCA), multi-dimensional scaling (MDS), Kruskal’s non-metric MDS (iMDS) (Kruskal, 1978), Sammon’s mapping (Sammon) (Sammon, 1969), locally linear embedding (LLE) (Roweis and Saul, 2000), Hessian LLE (HLLE) (Donoho and Grimes, 2003), isomap (Tenenbaum et al., 2000), kPCA, Laplacian eigenmap (LEIM), UMAP, t-SNE and PHATE (Moon et al., 2019). For methods such as kPCA, t-SNE, UMAP and PHATE, that require tuning parameters, we consider two different settings (Section A.2 of the Supplement) of tuning parameters for each method, denoted as “kPCA1” and “kPCA2,” etc. Therefore, for each dataset we obtain $K=16$ candidate visualizations corresponding to different combinations of visualization tools and tuning parameters. Applying our proposed method, we obtain eigenscores $\{\widehat{\mathbb{s}}_{i}\}_{1\leq i\leq n}$ for the candidate visualizations. We also compare two meta-distances based on the 16 visualizations, which are, the proposed spectral meta-distance matrix (“meta-spec”) based on the eigenscores, and the naive meta-distance matrix (“meta-aver”) assigning equal weights to all the candidate visualizations, as in (2.13).

To evaluate the proposed eigenscores, for each setting and each $i\in\{1,2,...,n\}$, we compute $\cos\angle(\widehat{\mathbb{s}}_{i},\mathbb{s}_{i}):=\frac{(\widehat{\mathbb{s}}_{i})^{\top}\mathbb{s}_{i}}{\|\widehat{\mathbb{s}}_{i}\|_{2}\|\mathbb{s}_{i}\|_{2}}$, for the angle between the eigenscores $\widehat{\mathbb{s}}_{i}$ and the true local concordance $\mathbb{s}_{i}$ defined as

where $\bar{\mathbb{P}}^{*}_{i.}$ is the $i$-th row of the normalized distance matrix $\bar{\mathbb{P}}^{*}$ for the underlying noiseless samples $\{\mathbb{Y}_{i}^{*}\}_{1\leq i\leq n}$, defined as in (2.5) with $\mathbb{X}_{i}^{(k)}$’s replaced by $\mathbb{Y}^{*}_{i}$’s. Table 1 shows empirical mean and standard error (SE) of the averaged cosines $\frac{1}{n}\sum_{i=1}^{n}\cos\angle(\widehat{\mathbb{s}}_{i},\mathbb{s}_{i})$, over the family of datasets under the same low-dimensional structure associated with various $\theta$ as shown in Figure 2(a). Our simulations showed that $\cos\angle(\widehat{\mathbb{s}}_{i},\mathbb{s}_{i})\approx 1$, indicating that the eigenscores $\widehat{\mathbb{s}}_{i}$ essentially characterize the true concordance between the patterns contained in each candidate visualization and that of the underlying noiseless samples, evaluated locally with respect to sample $i$. This justifies the proposed eigenscore as a precise measure of performance of the candidate visualizations in preserving the underlying true signals.

To assess the quality of two meta-distance matrices, for each dataset, we compare the mean concordance $\frac{1}{n}\sum_{i=1}^{n}(\bar{\mathbb{P}}_{i.}^{(k)})^{\top}\bar{\mathbb{P}}^{*}_{i.}$ between the normalized distance of each candidate visualization and that of the underlying noiseless samples, and the mean concordance $\frac{1}{n}\sum_{i=1}^{n}(\bar{\mathbb{P}}_{i.}^{m})^{\top}\bar{\mathbb{P}}^{*}_{i.}$ between the obtained meta-distance and that of the underlying noiseless samples. Figure 2(a) and Supplement Figure 7 show boxplots of these mean concordances for the 16 candidate visualizations and the two meta-distances under each setting of underlying structures across various values of $\theta$. We observe that for each of the three structures, our proposed meta-distance is substantially more concordant with the underlying true patterns, than every candidate visualization and the naive meta-distance, indicating the superiority of the proposed meta-distance. To further demonstrate the advantage of the spectral meta-distance and its benefits to the final meta-visualization, we compared our proposed meta-visualization using UMAP, and candidate visualizations of a dataset under setting (i) with $\theta=5$, and present their sample-wise concordance $\{(\bar{\mathbb{P}}_{i.}^{(k)})^{\top}\bar{\mathbb{P}}^{*}_{i.}\}_{1\leq i\leq n}$ for each $k$, and $\{(\bar{\mathbb{P}}_{i.}^{m})^{\top}\bar{\mathbb{P}}^{*}_{i.}\}_{1\leq i\leq n}$ of the proposed meta-distance (Figure 2(b) and Supplement Figure 8). We observe that, while each individual method may capture some clusters in the dataset but misses others, the proposed meta-visualization is able to combine strengths of all the candidate visualizations in order to capture all the underlying clusters. Finally, to demonstrate the flexibility of our method with respect to higher intrinsic dimension $r$, under the setting (i), we further evaluated the performance of different methods for $r\in\{15,30,50\}$. Supplement Figure 9 shows consistent and superior performance of the proposed method compared to the other approaches.

Cluster data are ubiquitous in scientific research and industrial applications. Our first real data example concerns $n=590$ fragments of text, extracted from English translations of eight religious books or sacred scripts including Book of Proverb (BOP), Book of Ecclesiastes (BOE1), Book of Ecclesiasticus (BOE2), Book of Wisdom (BOW), Four Noble Truth of Buddhism (BUD), Tao Te Ching (TTC), Yogasutras (YOG) and Upanishads (UPA) (Sah and Fokoué, 2019). All the text were pre-processed using natural language processing into a $590\times 8265$ Document Term Matrix that counts frequency of 8265 atomic words, such as truth, diligent, sense, power, in each text fragment. In other words, each text fragment was treated as a bag of words, represented by a vector with $8265$ features. The word counts were centred and normalized before downstream analysis.

As in our simulation studies, we still consider $K=16$ candidate visualizations generated by 12 different methods with various tuning parameters (see Section A.2 of the Supplement for implementation details). Figure 3(a) contains examples of candidate visualizations obtained by PHATE, t-SNE, and kPCA, whose median eigenscores were ranked top, middle and bottom among all the visualizations (Figure 3(b)), respectively. More examples are included in Supplement Figure 10. In each visualization, the samples (text fragments) were colored by their associated books, showing how well the visualization captures the underlying clusters of the samples. The usefulness and validity of the eigenscores in Figure 3(b) can be verified empirically, by visually comparing the clarity of cluster patterns demonstrated by each candidate visualizations in Figure 3(a) and in Supplement Figure 10. Figure 3(c) is the proposed meta-visualization¹¹1With a slight abuse of notation, we used “meta-spec” and “meta-aver” hereafter to refer to the final meta-visualizations rather than the meta-distance matrices as in Section 2.2.1. of the samples by applying UMAP to the meta-distance matrix, which shows substantially better clustering of the text fragments in accordance with their sources. In addition, the meta-visualization also reflected deeper relationship between the eight religious books, such as the similarity between the two Hinduism books YOG and UPA, the similarity between Buddhism (BUD) and Taoism (TTC), the similarity between the four Christian books BOE1, BOE2, BOP, and BOW, as well as the general discrepancy between Asian religions (Hinduism, Buddhism, Taoism) and non-Asian religions (Christianity). All of these important phenomena, while salient in our meta-visualization, only appeared vaguely in very few candidate visualizations such as those produced by PHATE (Figure 3(a)) and UMAP (Supplement Figure 10).

To quantitatively evaluate the preservation of the underlying clustering pattern, we computed for each visualization the Silhouette indices (Rousseeuw, 1987) with respect to the underlying true cluster membership, based on the normalized pairwise-distance matrices of the embeddings defined in (2.5). The Silhouette index (see Section A.2 of the Supplement for its definition), defined for each individual sample in a visualization, measures the amount of discrepancy between the within-class distances and the inter-class distances with respect to a given sample. As a result, for a given visualization, its Silhouette indices altogether indicate how well the underlying cluster pattern is preserved in a visualization, and higher Silhouette indices indicate that the underlying clusters are more separate. Empirically, we observed a notable correlation ($\rho=0.679$) between the median Silhouette indices and the median eigenscores across the candidate visualizations (Supplement Figure 11). In addition, for each candidate visualization, we found that samples with higher Silhouette index tend to have higher eigenscores (Supplement Figure 12), demonstrating the effectiveness of eigenscores, and its benefits on the final meta-visualization. In Figure 3(d), we show that, even taking into account the stochasticity of the visualization method (UMAP) applied to the meta-distance matrix, our meta-visualization had the median Silhouette index much higher than those of the candidate visualizations, as well as that of the meta-visualization “meta-aver” based on the naive meta-distance. It is of interest to note that “meta-spec” was the only visualization with a positive median Silhouette index, showing its better separation of clusters compared with other visualizations. Importantly, the proposed meta-visualization was not sensitive to the specific visualization method applied to the meta-distance matrices – similar results were obtained when we replaced UMAP by PHATE, the method having the highest median eigenscore in Figure 3(c), or t-SNE, for meta-visualization (Supplement Figure 11).

Our second real data example concerns visualization of a different low-dimensional structure, namely, a mixture of cycle and clusters, contained in the gene expression profile of a collection of mouse embryonic stem cells, as a result of the cell cycle mechanism. The cell cycle, or cell-division cycle, is the series of events that take place in a cell that cause it to divide into two daughter cells²²2https://en.wikipedia.org/wiki/Cell_cycle. Identifying the cell cycle stages of individual cells analyzed during development is important for understanding its wide-ranging effects on cellular physiology and gene expression profiles. Specifically, our dataset contains $n=288$ mouse embryonic stem cells, whose underlying cell cycle stages were determined using flow cytometry sorting (Buettner et al., 2015). Among them, one-third (96) of the cells are in the G1 stage, one-third in the S stage, and the rest in the G2M stage. The raw count data were preprocessed and normalized, leading to a dataset consisting of standardized expression levels of $1147$ cell-cycle related genes for the 288 cells (see Section A.2 of the Supplement for implementation details of our data preprocessing).

We obtained 16 candidate visualizations as before, and applied our proposed method. Figure 4(a) contains examples of candidate visualizations obtained by t-SNE, LEIM, and kPCA, whose median eigenscores were ranked top, middle and bottom among all the visualizations, respectively, and the cells were colored according to their true cell cycle stages. Figure 4(b) contains the boxplots of eigenscores for the candidate visualizations, indicating the overall quality of each visualization. The variation of eigenscores within each candidate visualization suggests that different visualizations have their own unique features and strengths to be contributed to the meta-visualization (Supplement Figure 16). Figure 4(c) is the proposed meta-visualization by applying kPCA to the meta-distance matrix. Comparing with Figure 4(a), the proposed meta-visualization showed better clustering of the cells according to their cell cycle stages, as well as a more salient cyclic structure underlying the three cell cycle stages (Supplement Figure 15). To quantify the performance of each visualization in terms of these two underlying structures (cluster and cycle), we considered two distinct metrics, namely, the median Silhouette index with respect to the underlying true cell cycle stages, and the Kendall’s tau statistic (Kendall, 1938) between the inferred relative order of the cells and their true orders on the cycle. Specifically, to infer the relative order of cells, we projected the coordinates of each visualization to the two-dimensional unit circle centred at the origin (Supplement Figure 15), and then determined the relative orders based on the cells’ respective projected positions on the unit circle. Figure 4(d) shows that the proposed meta-visualization was significantly better than all the candidate visualizations and the naive meta-visualization in representing both aspects of the data.

Our third real data example concerns visualization of a mixed pattern of a trajectory and clusters underlying the gene expression profiles of a collection of cells undergoing differentiation (Hayashi et al., 2018). Specifically, 421 mouse embryonic stem cells were induced to differentiate into primitive endoderm cells. After the induction of differentiation, the cells were dissociated and individually captured at 12- or 24-hour intervals (0, 12, 24, 48 and 72 h), and each cell was sequenced to obtain the final total RNA sequencing reads using the random displacement amplification sequencing technology. As a result, at each of the five time points, there were about 70 to 90 cells captured and sequenced. The raw count data were preprocessed and normalized (see Section A.2 of the Supplement for implementation details), leading to a dataset consisting of standardized expression levels of $500$ most variable genes for the 421 cells.

Again, we obtained 16 candidate visualizations as before, and applied our proposed method. In Figure 5(a)-(c) we show examples of candidate visualizations, boxplots of the eigenscores, and the meta-visualization using kPCA. The global (Figure 5(b)) and local (Supplement Figure 18) variation of eigenscores demonstrated contribution of different visualizations to the final meta-visualization according to their respective performance. We observed that some candidate visualizations such as kPCA, UMAP (Figure 5(a)) and PHATE (Supplement Figure 17) to some extent captured the underlying trajectory structure consistent with the time course of the cells. However, the meta-visualization in Figure 5(c) showed much more salient patterns in terms of both the underlying trajectory and the cluster pattern among the cells, by locally combining strengths of the individual visualizations (Supplement Figure 18). We quantified the performance of visualizations from these two aspects using the median Silhouette index with respect to the underlying true cluster membership (i.e., batches of time course) and Kendall’s tau statistic between the inferred cell order and the true order along the progression path. To infer the relative order of the cells from a visualization, we ordered all the cells based on the two-dimensional embedding along the direction that explained the most variability of the cells. In Figure 5(d), we observed that, the proposed meta-visualization had the largest median Silhouette index as well as the largest Kendall’s tau statistic, compared with all the candidate visualizations and the naive meta-visualization, showing the superiority of the proposed meta-visualization in both aspects.

For datasets of moderate size as the ones analyzed in the previous sections, the proposed method had a computational cost comparable to that of t-SNE or UMAP for generating a single candidate visualization (Supplement Figure 19). As for very large and high-dimensional datasets, there are a few features of the proposed algorithm that make it readily scalable. First, although our method relies on computing the leading eigenvector of generally non-sparse matrices, these matrices (i.e., $\mathbb{G}_{i}$ in Algorithm 1) are of dimension $K\times K$, where $K$ – the number of candidate visualizations – is usually much smaller compared to the sample size $n$ or dimensionality $p$ of the original data. Thus, for each sample $i$, the computational cost due to the eigendecomposition is mild. Second, given the candidate visualizations, our proposed algorithm is independent of the dimensionality ($p$) of the original dataset, as it only requires as input a set of low-dimensional embeddings produced by different visualization methods. Third, since our algorithm computes the eigenscores and the meta-distance with respect to each sample individually, the algorithm can be easily parallelized and carried out in multiple cores to further reduce time cost.

To demonstrate the computational efficiency of the proposed method for large and high-dimensional datasets, we evaluated the proposed method on real single-cell transcriptomic datasets (Buckley et al., 2022) of various sample sizes ($n\in\{1000,2000,4000,8000,14000\}$ cells of nine different cell types from the neurogenic regions of mice) and dimensions ($p\in\{500,1000,2000\}$ genes). For each dataset, we obtained 11 candidate visualizations and applied Algorithm 1 to generate the final meta-visualization (see Section A.2 of the Supplement for implementation details). Supplement Figure 20(b) contains boxplots of median Silhouette indices for each candidate visualizations and the meta-visualization (highlighted in red) with respect to the underlying true cell types, showing the stable and superior performance of the proposed method under various sample sizes and dimensions. In Supplement Figure 20(a), we compared the running time for generating the 11 candidate visualizations, and that for generating the meta-visualizations based on Algorithm 1, on a MacBook Pro with 2.2 GHz 6-Core Intel Core i7. In general, as $n$ became large, the running time of the proposed algorithm also increased, but remained much less than that for generating the candidate visualizations. The difference in time cost became more significant as $n$ increased, demonstrating that for very large and high-dimensional datasets the computational cost essentially comes from generating candidate visualizations, rather than from the meta-visualization step. In particular, for dataset of sample size as large as $8000$ and of dimension $2000$, it took about 60 mins to generate all the 11 candidate visualizations, and took about additional 12 mins to generate the meta-visualization. Moreover, Supplement Figure 20(a) also demonstrated that, for each $n$, when $p$ increased, the running time for generating the candidate visualizations was longer, but the time cost for meta-visualization remained about the same (difference less than one minute). We also note that users often create multiple visualizations for data exploration, and our approach can simply reuse these visualizations with little additional computational cost.

We develop a general and flexible theoretical framework, to investigate the statistical properties of the proposed methods, as well as the fundamental principles behind its empirical success³³3Throughout, we adopt the following notations. For a matrix $\mathbb{A}=(a_{ij})\in\mathbb{R}^{n\times n}$, we define its spectral norm as $\|\mathbb{A}\|=\sup_{\|\mathbb{x}\|_{2}\leq 1}\|\mathbb{A}\mathbb{x}\|_{2}$. For sequences $\{a_{n}\}$ and $\{b_{n}\}$, we write $a_{n}=o(b_{n})$ or $b_{n}\gg a_{n}$ if $\lim_{n}a_{n}/b_{n}=0$, and write $a_{n}\asymp b_{n}$ if there exists constants $C_{1},C_{2}>0$ such that $C_{1}b_{n}\leq a_{n}\leq C_{2}b_{n}$ for all $n$.. As can be seen from Section 2.1, there are two key ingredients of our proposed method, namely, the eigenscores $\{\widehat{\mathbb{s}}_{i}\}_{1\leq i\leq n}$ for evaluating the candidate visualizations, and the meta-distance matrix $\bar{\mathbb{P}}^{m}$ that combines multiple candidate visualizations to obtain a meta-visualization.

To formally study their properties, we introduce a generic model for the collection of $K$ candidate visualizations produced by multiple visualization methods, with possibly different settings of tuning parameters for a single method as considered in previous sections. Specifically, we assume $\{\mathbb{Y}_{i}\}_{1\leq i\leq n}$ are generated as

where $\{\mathbb{Y}_{i}^{*}\}_{1\leq i\leq n}$ are the underlying noiseless samples and $\{\mathbb{Z}_{i}\}_{1\leq i\leq n}$ are the random noises. Recall that $\{\mathbb{P}^{(k)}\}_{1\leq k\leq K}$ are the distance matrices associated to the candidate visualizations. Then, for the candidate visualizations, we consider a scaled signal-plus-noise expression

induced by (2.11), where $c_{i,k}\geq 0$ is a global scaling parameter, $\mathbb{P}^{*}_{i.}$ is the $i$-th row of the pairwise distance matrix $\mathbb{P}^{*}=(\|\mathbb{Y}^{*}_{i}-\mathbb{Y}^{*}_{j}\|_{2})_{1\leq i,j\leq n}$ of the underlying noiseless samples, and $\mathbb{h}_{i}^{(k)}$ is a random vector characterizing the relative distortion of $\mathbb{P}_{i.}^{(k)}$ associated to the $k$-th candidate visualization, from the underlying true pattern $\mathbb{P}^{*}_{i.}$. Before characterizing the distributions of $\{\mathbb{h}_{i}^{(k)}\}_{1\leq k\leq K}$, we point out that, in principle, the relative distortions $\{\mathbb{h}_{i}^{(k)}\}_{1\leq k\leq K}$ are jointly determined by the random noises $\{\mathbb{Z}_{i}\}_{1\leq i\leq n}$ in (2.11), and the features and relations between of the specific visualization methods. Importantly, in line with what is often encountered in practice, equation (2.12) allows for flexible and possibly distinct scaling and directionality for different candidate visualizations, by introducing the visualization-specific parameter $c_{k}$, and by focusing on the pairwise distance matrices, rather than the low-dimensional embeddings $\{\mathbb{X}_{i}^{(k)}\}_{1\leq i\leq n}$ themselves.

To quantitatively describe the variability of the distortions $\{\mathbb{h}_{i}^{(k)}\}_{1\leq k\leq K}$ across $K$ candidate visualizations, we assume

• (C1a)

  $\{\mathbb{h}_{i}^{(k)}\}_{1\leq k\leq K}$ are identically distributed sub-Gaussian vectors with parameter $\sigma^{2}$, that is, for any deterministic unit vector $\mathbb{g}\in\mathbb{R}^{n}$, we have $\mathbb{E}\exp\{(\mathbb{h}^{(k)}_{i})^{\top}\mathbb{g}\}\leq\exp(\sigma^{2}/2)$, and that $\|\mathbb{h}^{(k)}_{i}\|^{2}_{2}=c\sigma^{2}n(1+o(1))$ with high probability⁴⁴4An event $A_{n}$ holds with high probability if there exists some $N>0$, such that $P(A_{n})\geq 1-n^{-D}$ for all $n\geq N$ for some large constant $D>0$. for some constant $c>0$.

This assumption makes (2.12) a generative model for $\{\mathbb{P}_{i.}^{(k)}\}_{1\leq k\leq K}$ with ground truth $\mathbb{P}^{*}_{i.}$ and random distortions, where the variance parameter $\sigma$ describes the average level of the distortions of candidate visualizations from the truth after proper scaling. In relation to (2.11), such a condition can be satisfied when the signal structure $\{\mathbb{Y}_{i}^{*}\}_{1\leq i\leq n}$ is finite, the noise $\{\mathbb{Z}_{i}\}_{1\leq i\leq n}$ is sub-Gaussian, and the dimension reduction map underlying the candidate visualization is bounded and sufficiently smooth. See Section B.2 of the Supplement for details. In addition, we also need to characterize the correlations among these random distortions, not only because the candidate visualizations are typically obtained from the same dataset $\{\mathbb{Y}_{i}\}_{1\leq i\leq n}$, but also because of the possible similarity between the adopted visualization methods, such as MDS and iMDS, or t-SNE under different tuning parameters. Specifically, for any $j,k\in\{1,2,...,K\}$, we define the cross-visualization covariance $\mathbb{\Sigma}_{jk}=\mathbb{E}\mathbb{h}^{(j)}_{i}(\mathbb{h}^{(k)}_{i})^{\top}$, and quantify the level of dependence between a pair of candidate visualizations by $\rho_{jk}=\|\mathbb{\Sigma}_{jk}\|/\sigma^{2}$. By Condition (C1a), we have $\rho_{jj}\leq 1$ for all $j$. For all correlation parameters $\{\rho_{jk}\}_{1\leq j,k\leq K}$, we assume

• (C1b)

  The matrix $\mathbb{R}=(\rho_{jk})_{1\leq j,k\leq K}$ satisfies $\rho:=\|\mathbb{R}\|=o(K)$.

Condition (C1b) covers a wide range of correlation structures among the candidate visualizations, allowing in particular for a subset of highly correlated visualizations possibly produced by very similar methods. The parameter $\rho$ characterizes the overall correlation strength among the candidate visualizations, which is assumed to be not too large. As a comparison, note that a set of pairwise independent candidate visualizations implies that $\rho\approx 1$, whereas a set of identical candidate visualizations have $\rho\approx K$. In particular, the requirement $\rho=o(K)$ can be satisfied if, for example, among $K$ candidate visualizations, there are subsets of at most $\sqrt{K}$ visualizations that are produced by very similar procedures, such as by the same method under different tuning parameters, so that $\rho\leq\sqrt{K}=o(K)$. When Condition (C1b) fails, as all the candidate visualizations are essentially similarly distorted from truth, combination of them will not be substantially more informative than each individual visualization.

Under Condition (C1a), it holds that $\mathbb{E}\|\mathbb{h}_{i}^{(k)}\|_{2}\asymp\sigma\sqrt{n}$. Hence, we can use the quantity $\frac{\|\mathbb{P}^{*}_{i.}\|_{2}}{\sigma\sqrt{n}}$ to characterize the overall SNR in the candidate visualizations as modelled by (2.12), which reflects the average quality of the candidate visualizations in preserving the underlying true patterns around sample $i$. Before stating our main theorems, we first introduce our main assumption on the minimal SNR requirement, that is,

• (C2)

  For $(\sigma,\rho)$ defined in (C1a) and (C2b), it holds that $\frac{\|\mathbb{P}^{*}_{i.}\|_{2}}{\sigma\sqrt{n}}\gg\sqrt{\rho/K}$ and $K=o(n)$ as $n\to\infty$.

Our algorithm is expected to perform well if $\sqrt{\rho/K}$ is small relative to the overall SNR. The condition $K=o(n)$ is easily satisfied for a sufficiently large dataset.

Recall that $\mathbb{s}_{i}=((\bar{\mathbb{P}}_{i.}^{(1)})^{\top}\bar{\mathbb{P}}^{*}_{i.},(\bar{\mathbb{P}}_{i.}^{(2)})^{\top}\bar{\mathbb{P}}^{*}_{i.},...,(\bar{\mathbb{P}}_{i.}^{(K)})^{\top}\bar{\mathbb{P}}^{*}_{i.}).$ The following theorem concerns the convergence of eigenscores to the true concordance $\mathbb{s}_{i}$, and is proved in Section B.3 of the Supplement.

Theorem 1 implies that, as long as the candidate visualizations contain sufficient amount of information about the underlying true structure, and are not terribly correlated, the proposed eigenscores $\{\widehat{\mathbb{s}}_{i}\}_{1\leq i\leq n}$ are quantitatively reliable, as they converge to the actual quality measures $\{\mathbb{s}_{i}\}_{1\leq i\leq n}$ asymptotically. In other words, the eigenscores provide a point-wise consistent estimation of the concordance between the candidate visualizations as summarized by $\{\mathbb{P}^{(k)}\}_{1\leq k\leq K}$ and the underlying true patterns $\mathbb{P}^{*}$, justifying the empirical observations in Table 1. Importantly, Condition (C2) suggests that our proposed eigenscores may benefit from a larger number $K$ of candidate visualizations, or a smaller overall correlation $\rho$, that is, a collection of functionally more diverse candidate visualizations; see further discussions in Section 2.3.6.

Our second theorem concerns the guaranteed performance of our proposed meta-distance matrix and its improvement upon the individual candidate visualizations in the large-sample limit.

Theorem 2 is proved in Section B.4 of the Supplement. In addition to the point-wise consistency of $\bar{\mathbb{P}}^{m}$ as described by $\cos\angle(\bar{\mathbb{P}}^{m}_{i.},\mathbb{P}^{*}_{i.})\to 1$ in probability, Theorem 2 also ensures that the proposed meta-distance is in general no worse than the individual candidate visualizations, suggesting a competitive performance of the meta-visualization. In particular, if in addition to Conditions (C1a) (C1b) and (C2) we also have $\frac{\|\mathbb{P}_{i.}^{*}\|_{2}}{\sigma\sqrt{n}}\leq C$, that is, the magnitude of the random distortions from the true structure $\mathbb{P}_{i.}^{*}$ is relatively large, then each candidate visualization necessarily has at most mediocre performance, i.e., $\max_{1\leq k\leq K}\cos\angle(\mathbb{P}^{(k)}_{i.},\mathbb{P}^{*}_{i.})<1-\delta$ in probability. In such cases, the proposed meta-distances is still consistent and thus strictly better than all candidate visualizations. Theorem 2 justifies the superior performance of the spectral meta-visualization demonstrated in Section 2.2, compared with 16 candidate visualizations.

Among the three conditions required for the consistency of the proposed meta-distance matrix, Condition (C2) is most critical as it describes the minimal SNR requirement, that is, how much information the candidate visualizations altogether should contain about the underlying true structure of the data. In this connection, our theoretical analysis indicates that, in fact, such a signal strength condition is also necessary, not only for the proposed method, but for any possible methods. More specifically, in Section B.6 of the Supplement, we proved (Theorem 4) that, it’s impossible to construct a meta-distance matrix that is consistent when Condition (C2) is violated. This result shows that the settings where our meta-visualization algorithm works well is essentially the most general setting possible.

In Section 2.2, in addition to the proposed meta-visualization, we also considered the meta-visualization based on the naive meta-distance matrix $\bar{\mathbb{P}}^{a}$, whose rows are

which is a simple average across all the candidate visualizations. We observed in all our real-world data analyses that, such a naive meta-visualization only had mediocre performance compared to the candidate visualizations (Figures 3, 4 and 5), much worse than the proposed spectral meta-visualization. The empirical observations suggest the advantage of informative weighting for combining candidate visualizations.

The empirically observed suboptimality of the non-informative weighting procedure can justified rigorously by theory. Our next theorem concerns the behavior of the proposed meta-distance matrix $\bar{\mathbb{P}}^{m}$ and the naive meta-distance matrix $\bar{\mathbb{P}}^{a}$ when combining a mixture of well-conditioned candidate visualizations, as characterized by our assumptions (C1a) (C1b) and (C2), and some adversarial candidate visualizations whose pairwise-distance matrices does not contain any information about the true structure. Specifically, we suppose among all the $K$ candidate visualizations, there is a collection $\mathcal{C}_{0}$ of $(1-\eta)K$ well-conditioned candidate visualizations for some small $\eta\in(0,1)$, and a collection $\mathcal{C}_{1}$ of $\eta K$ adversarial candidate visualizations.

Theorem 3 is proved in Section B.5 of the Supplement. By Theorem 3, on the one hand, even when there are a small portion of really poor (adversarial) candidate visualizations to be combined with other relatively good visualizations, the proposed method still perform well thanks to the consistent eigenscore weighting in light of Theorem 1. On the other hand, no matter how strong the SNR is for those well-conditioned candidate visualizations, the method based on non-informative weighting is strictly sub-optimal. Indeed, when $\|\mathbb{P}_{i.}^{*}\|_{2}\gg\sigma\sqrt{n}\asymp\mathbb{E}\|\mathbb{h}_{i}\|_{2}$, although we have $\cos\angle(\mathbb{P}_{i.}^{(k)},\mathbb{P}^{*}_{i.})\to 1$ in probability for all $k\in\mathcal{C}_{0}$, the non-informative weighting would suffer from the non-negligible negative effects from the adversarial visualizations in $\mathcal{C}_{1}$, causing a strict deviation from $\mathbb{P}^{*}$; see, for example, Figures 3-5(b)(d) for empirical evidences from real-world data.

The proposed eigenscores provide an efficient and consistent way of evaluating the performance of the candidate visualizations. As mentioned in Section 1, a number of metrics have been proposed to quantify the distortion of a visualization by comparing the low-dimensional embedding directly with the original high-dimensional data. Such metrics essentially treat the original high-dimensional data as the ground truth, and do not take into account the noisiness of the high-dimensional data. However, for many datasets arising from real-world applications, the observed datasets, as modelled by (2.11), are themselves very noisy, which may not make an ideal reference point for evaluating a visualization that probably has already significantly denoised the data through dimension reduction. For example, the three real-world datasets considered in Section 2.2 are all high-dimensional and contain much more features than number of samples. In each case, there are some underlying clusters among the samples, but the original datasets showed significantly weaker cluster structure compared to most of the 16 candidate visualizations (Supplement Figure 21), suggesting that directly comparing a visualization with the noisy high-dimensional data may be misleading. In this respect, our theorems indicate that the proposed spectral method is able to precisely assess and effectively combine multiple visualizations to better grasp the underlying noiseless structure $\mathbb{P}^{*}$, without referring to the original noisy datasets, making it more robust, flexible, and computationally more efficient.

Our theoretical analysis implies that the proposed meta-visualization may benefit from a large number (larger $K$) of functionally diverse (small $\rho$) candidate visualizations. To empirically verify this theoretical observation, we focused on the religious and biblical text data and the mouse embryonic stem cells data considered in Section 2.2, and obtained spectral meta-visualizations based on a smaller but relatively diverse collection of 5 candidate visualizations, produced by arguably the most popular methods, namely, t-SNE, PHATE, UMAP, PCA and MDS, respectively. Compared with the 16 candidate visualizations considered in Section 2.2, here we have presumably similar $\rho$ but much smaller $K$. As a result, for the religious and biblical texts data, the meta-visualization had a median Silhouette index 0.187 (Supplement Figure 13), which was smaller than the median Silhouette index 0.275 based on the 16 candidate visualizations as in Figure 3(d); for the cell cycle data, the meta-visualization had a median Silhouette index -0.062 and a Kendall’s tau statistic 0.313, both smaller than the respective values based on the 16 candidate visualizations as in Figure 4(d). On the other hand, we also evaluated the effect when $\rho$ is increased but $K$ remains fixed. Specifically, we obtained 16 candidate visualizations, all produced by PHATE with varying nearest neighbor parameters, the final spectral meta-visualization had a median Silhouette index 0.094, which was even lower than the above meta-visualization based on five distinct methods, although being still slightly better than the 16 PHATE-based candidate visualizations (Supplement Figure 13). These empirical evidences were in line with our theoretical predictions, suggesting benefits of including more diverse visualizations.