Multi-view data visualisation via manifold learning

Notation

Throughout this article, the following notation is used:

• N: The number of samples.

• $p$: The number of variables of the design matrix.

• $X\in {\mathbb{R}}^{N\times p}$: A single-view data matrix, representing the original high-dimensional data used as input; ${\mathbf{x}}_i\in {\mathbb{R}}^p$ is the $i^{th}$ data point of X.

• M: The number of data-views in a given data set; $m\in \{1,\cdots ,M\}$ represents an arbitrary data-view.

• $X^{(m)}\in {\mathbb{R}}^{N\times p_m}$: The $m^{th}$ data-view of multi-view data; ${\mathbf{x}}_i^m\in {\mathbb{R}}^{p_m}$ is the $i^{th}$ data point of $X^{(m)}$.

• $Y\in {\mathbb{R}}^{N\times d}$: A low-dimensional embedding of the original data. ${\mathbf{y}}_i\in {\mathbb{R}}^d$ represents the $i^{th}$ data point of Y. In this manuscript, $d=2$, as the focus of the manuscript is on data visualisation.

Multi-SNE

SNE, proposed by Hinton & Roweis (2003), measures the probability distribution, P of each data point ${\mathbf{x}}_i$ by looking at the similarities among its neighbours. For every sample $i$ in the data, $j$ is taken as its potential neighbour with probability $p_{ij}$, given by

(1) $$p_{ij}=\frac{\exp (-d_{ij}^2)}{{\sum }_{k\ne i}\exp (-d_{ik}^2)},$$where $d_{ij}=\frac{||{\mathbf{x}}_i-{\mathbf{x}}_j|{|}^2}{2{\sigma }_i^2}$ represents the dissimilarity between points ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$. The value of ${\sigma }_i$ is either set by hand or found by binary search (van der Maaten & Hinton, 2008). Based on this value, a probability distribution of sample $i$, $P_i={\sum }_j{p}_{ij}$, with fixed perplexity is produced. Perplexity refers to the effective number of local neighbours and it is defined as $Perp(P_i)=2^{H(P_i)}$, where $H(P_i)=-{\sum }_j{p}_{ij}{\log }_2{p}_{ij}$ is the Shannon entropy of $P_i$. It increases monotonically with the variance ${\sigma }_i$ and typically takes values between $5$ and $50$.

In the same way, a probability distribution in the low-dimensional space, Y, is computed as follows:

(2) $$q_{ij}=\frac{\exp (-||{\mathbf{y}}_i-{\mathbf{y}}_j|{|}^2)}{{\sum }_{k\ne i}\exp (-||{\mathbf{y}}_i-{\mathbf{y}}_k|{|}^2)},$$which represents the probability of point $i$ selecting point $j$ as its neighbour.

The induced embedding output, ${\mathbf{y}}_i$, represented by probability distribution, Q, is obtained by minimising the Kullback-Leibler divergence (KL-divergence) $KL(P||Q)$ between the two distributions P and Q (Kullback & Leibler, 1951). The aim is to minimise the cost function:

(3) $$C_{SNE}={\sum }_{i}KL(P_i||Q_i)={\sum }_i{\sum }_j{p}_{ij}\log \frac{p_{ij}}{q_{ij}}$$

Hinton & Roweis (2003) assumed a Gaussian distribution in computing the similarity between two points in both high and low dimensional spaces. van der Maaten & Hinton (2008) proposed a variant of SNE, called t-SNE, which uses a symmetric version of SNE and a Student t-distribution to compute the similarity between two points in the low-dimensional space Q, given by

(4) $$q_{ij}=\frac{{(1+||{\mathbf{y}}_i-{\mathbf{y}}_j|{|}^2)}^{-1}}{{\sum }_{k\ne l}{(1+||{\mathbf{y}}_k-{\mathbf{y}}_l|{|}^2)}^{-1}}$$

T-SNE is often preferred, because it reduces the effect of crowding problem (limited area to accommodate all data points and differentiate clusters) and it is easier to optimise, as it provides simpler gradients than SNE (van der Maaten & Hinton, 2008).

We propose multi-SNE, a multi-view manifold learning algorithm based on t-SNE. Our proposal computes the KL-divergence between the distribution of a single low-dimensional embedding and each data-view of the data separately, and minimises their weighted sum. An iterative algorithm is proposed, in which at each iteration the induced embedding is updated by minimising the cost function:

(5) $$C_{multi-SNE}={\sum }_m{\sum }_i{\sum }_j{w}^m{p}_{ij}^m\log \frac{p_{ij}^m}{q_{ij}},$$where $w^m$ is the combination coefficient of the $m^{th}$ data-view. The vector $\mathbf{w}=(w^1,\cdots ,w^M)$ acts as a weight vector that satisfies ${\sum }_m{w}^m=1$. In this study, equal weights on all data-views were considered, i.e. $w^m=\frac{1}{M},\mathrm{\forall }m=1,\cdots ,M$. The algorithm of the proposed multi-SNE approach is presented in Algorithm 1 of Section 1 in the Supplemental File.

An alternative multi-view extension of t-SNE, called m-SNE was proposed by Xie et al. (2011). M-SNE applies t-SNE on a single distribution in the high-dimensional space, which is computed by combining the probability distributions of the data-views, given by $p_{ij}={\sum }_{m=1}^M{\beta }^m{p}_{ij}^m$. The coefficients (or weights) ${\beta }^m$ share the same role as $w^m$ in multi-SNE and similarly $\text{β}=({\beta }^1,\cdots ,{\beta }^M)$ satisfies ${\sum }_m{\beta }^m=1$. This leads to a different cost function than the one in Eq. (5).

Kanaan Izquierdo (2017) proposed a similar cost function for multi-view t-SNE, named MV-tSNE1, given as follows:

(6) $$C_{\mathrm{M}\mathrm{V}-\mathrm{t}\mathrm{S}\mathrm{N}\mathrm{E}1}=\sum _m\sum _i\sum _j{p}_{i|j}^m\log \frac{p_{i|j}^m}{q_{i|j}}$$

Their proposal is a special case of multi-SNE, with $w_m=\frac{1}{M}$. Kanaan Izquierdo (2017) did not pursue MV-tSNE1 any further, but instead, they proceeded with an alternative solution, MV-tSNE2, which combines the probability distributions (similar to m-SNE) through expert opinion pooling. A comparison between multi-SNE, m-SNE and MV-tSNE2 is presented in Fig. S2 of Section 4.1 in the Supplemental File. Based on two real data sets, multi-SNE and m-SNE outperformed MV-tSNE2, with the solution by multi-SNE producing the best separation among the clusters in both examples.

Multi-SNE avoids combining the probability distributions of all data-views together. Instead, the induced embeddings are updated by minimising the KL-divergence between every data-view’s probability distribution and that of the low-dimensional representation we seek to obtain. In other words, this is achieved by computing and summing together the gradient descent for each data-view. The induced embedding is then updated by minimising the summed gradient descent.

Throughout this article, for all variations of t-SNE we have applied the PCA pre-training step proposed by van der Maaten & Hinton (2008). van der Maaten & Hinton (2008) discussed that by reducing the dimensions of the input data through PCA the computational time of t-SNE is reduced. In this article, the principal components taken retained at least $80\mathrm{\%}$ of the total variation (variance explained) in the original data. In addition, as the multi-SNE algorithm is an iterative algorithm we opted for running the algorithm for 1,000 iterations for all analyses conducted. Alternatively, a stopping rule could have been implemented with the iterative algorithm to stop after no significant changes were observed to the cost-function. Both these options are available at the implementation of the multi-SNE algorithm.

Multi-LLE

LLE attempts to discover a non-linear structure of high-dimensional data, X, by computing low-dimensional and neighbourhood-preserving embeddings, Y (Saul & Roweis, 2001). The main three steps of the algorithm are:

• 1.

  The set, denoted by ${\mathrm{\Gamma }}_i$, contains the K nearest neighbours of each data point ${\mathbf{x}}_i,i=1,\cdots ,N$. The most common distance measure between the data points is the Euclidean distance. Other local metrics can also be used in identifying the nearest neighbours (Roweis & Saul, 2000).

• 2.

  A weight matrix, W, is computed, which acts as a bridge between the high-dimensional space in X and the low-dimensional space in Y. Initially, W reconstructs X, by minimising the cost function:

  (7) $${\mathcal{E}}_X=\sum _i|{\mathbf{x}}_i-\sum _j{W}_{ij}{\mathbf{x}}_j{|}^2$$where the weights $W_{ij}$ describe the contribution of the $j^{th}$ data point to the $i^{th}$ reconstruction. The optimal weights $W_{ij}$ are found by solving the least squares problem given in Eq. (7) subject to the constraints:

  • (a) $W_{ij}=0$, if $j\notin {\mathrm{\Gamma }}_i$, and

  • (b) ${\sum }_j{W}_{ij}=1$

• 3.

  Once W is computed, the low-dimensional embedding ${\mathbf{y}}_i$ of each data point $i=1,\cdots ,N$, is obtained by minimising:

  (8) $${\mathcal{E}}_Y=\sum _i|{\mathbf{y}}_i-\sum _j{W}_{ij}{\mathbf{y}}_j{|}^2$$

  The solution to Eq. (8), is obtained by taking the bottom $d$ non-zero eigenvectors of the sparse $N\times N$ matrix, $M=(I-W{)}^T(I-W)$ (Roweis & Saul, 2000).

We propose multi-LLE, a multi-view extension of LLE, that computes the low-dimensional embeddings by using the consensus weight matrix:

(9) $$\hat{W}=\sum _m{\alpha }^m{W}^m$$where ${\sum }_m{\alpha }^m=1$, and $W^m$ is the weight matrix for each data-view $m=1,\cdots M$. Thus, $\hat{Y}$ is obtained by solving:

$${\mathcal{E}}_{\hat{Y}}=\sum _i|{\hat{\mathbf{y}}}_i-\sum _j{\hat{W}}_{ij}{\hat{\mathbf{y}}}_j{|}^2$$

The multi-LLE algorithm is presented in Algorithm 2 of Section 1 in the Supplemental File.

Shen, Tao & Ma (2013) proposed m-LLE, an alternative multi-view extension of LLE. The LLE embeddings of each data-view are combined and LLE is applied to each data-view separately. The weighted average of those embeddings is taken as the unified low-dimensional embedding. In other words, computing the weight matrices $W^m$ and solving ${\mathcal{E}}_{Y^m}={\sum }_i|{\mathbf{y}}_i^m-{\sum }_j{W}_{ij}^m{\mathbf{y}}_j^m{|}^2$, for each $m=1,\cdots M$ separately. Thus, the low-dimensional embedding $\hat{Y}$ is computed by $\hat{Y}={\sum }_m{\beta }^m{Y}^m$, where ${\sum }_m{\beta }^m=1$.

An alternative multi-view LLE solution was proposed by Zong et al. (2017) to find a consensus manifold, which is then used for multi-view clustering via non-negative matrix factorization; we refer to this approach as MV-LLE. This solution minimises the cost function by assuming a consensus weight matrix across all data-views, as given in Eq. (9). The optimisation is then solved by using the Entropic Mirror Descent Algorithm (EMDA) (Beck & Teboulle, 2003). In contrast to m-LLE and MV-LLE, multi-LLE combines the weight matrices obtained from each data-view, instead of the LLE embeddings. No comparisons were conducted between MV-LLE and the proposed multi-LLE, as the code of the MV-LLE algorithm is not publicly available.

Multi-ISOMAP

ISOMAP aims to discover a low-dimensional embedding of high-dimensional data by maintaining the geodesic distances between all points (Tenenbaum, Silva & Langford, 2000); it is often regarded as an extension of Multi-dimensional Scaling (MDS) (Kruskal, 1964). The ISOMAP algorithm comprises of the following three steps:

• Step 1.

  A graph is defined. Let $G\sim (V,E)$ define a neighbourhood graph, with vertices V representing all data points. The edge length between any two vertices $i,j\in V$ is defined by the distance metric $d_X(i,j)$, measured by the Euclidean distance. If a vertex $j$ does not belong to the K nearest neighbours of $i$, then $d_X(i,j)=\mathrm{\infty }$. The parameter K is given as input, and it represents the connectedness of the graph G; as K increases, more vertices are connected.

• Step 2.

  The shortest paths between all pairs of points in G are computed. The shortest path between vertices $i,j\in V$ is defined by $d_G(i,j)$. Let $D_G\in {\mathbb{R}}^{|V|\times |V|}$ be a matrix containing the shortest paths between any vertices $i,j\in V$, defined by $(D_G{)}_{ij}=d_G(i,j)$.

  The most efficient known algorithm to perform this task is Dijkstra’s Algorithm (Dijkstra, 1959). In large graphs, an alternative approach to Dijkstra’s Algorithm would be to initialize $d_G(i,j)=d_X(i,j)$ and replace all entries by $d_G(i,j)=min\left\{d_G(i,k),d_G(k,j)\right\}$.

• Step 3.

  The low-dimensional embeddings are constructed. The $i^{th}$ component of the low-dimensional embedding is given by $y_i=\sqrt{{\lambda }_p}{u}_p^i$, where $u_p^i$ the $i^{th}$ component of $p^{th}$ eigenvector and ${\lambda }_p$ is the $p^{th}$ eigenvalue in decreasing order of the the matrix $\tau (D_G)$ (Tenenbaum, Silva & Langford, 2000). The operator, $\tau$ is defined by $\tau (D)=-\frac{HSH}{2}$, where S is the matrix of squared distances defined by $S_{ij}=D_{ij}^2$, and H is defined by $H_{ij}={\delta }_{ij}-\frac{1}{N}$. This is equivalent to applying classical MDS to $D_G$, leading to a low-dimensional embedding that best preserves the manifold’s estimated intrinsic geometry.

Multi-ISOMAP is our proposal for adapting ISOMAP on multi-view data. Let $G_m\sim (V,E_m)$ be a neighbourhood graph obtained from data-view $X^{(m)}$ as defined in the first step of ISOMAP. All neighbourhood graphs are then combined into a single graph, $\tilde{G}$; this combination is achieved by computing the edge length as the averaged distance of each data-view, i.e. $d_{\tilde{G}}(i,j)=w_m{\sum }_m{d}_{G_m}(i,j)$. Once a combined neighbourhood graph is computed, multi-ISOMAP follows steps 2 and 3 of ISOMAP described above. For simplicity, the weights throughout this article were set as $w_m=\frac{1}{M},\mathrm{\forall }m$. The multi-ISOMAP algorithm is presented in Algorithm 3 of Section 1 in the Supplemental File.

For completion, we have in addition adapted ISOMAP for multi-view visualisation following the framework of both m-SNE and m-LLE. Following the same logic, m-ISOMAP combines the ISOMAP embeddings of each data-view by taking the weighted average of those embeddings as the unified low-dimensional embedding. In other words, the low-dimensional embedding $\hat{Y}$ is obtained by computing $\hat{Y}={\sum }_m{\beta }^m{Y}^m$, where ${\sum }_m{\beta }^m=1$.

Parameter tuning

The multi-view manifold learning algorithms were tested on real and synthetic data sets for which the samples can be separated into several clusters. The true clusters are known and they were used to tune the parameters of the methods. To quantify the clustering performance, we used the following four extrinsic measures: (i) accuracy (ACC), (ii) Normalised Mutual Information (NMI) (Vinh, Epps & Bailey, 2010), (iii) Rand Index (RI) (Rand, 1971) and (iv) Adjusted Rand Index (ARI) (Hubert & Arabie, 1985). All measures take values in the range $\left[0,1\right]$, with $0$ expressing complete randomness, and $1$ perfect separation between clusters. The mathematical formulas of the four measures are presented in Section 2 of the Supplemental File.

SNE, LLE and ISOMAP depend on parameters of which their proper tuning ensues to optimal results. LLE and ISOMAP depend on the number of nearest neighbours (NN). SNE depends on the Perplexity ( $Perp$) parameter, which is directly related to the number of nearest neighbours. Similarly, the multi-view extensions of the three methods depend on the same parameters. The choice of the parameter can influence the visualisations and in some cases present the data into separate maps (van der Maaten & Hinton, 2012).

By assuming that the data samples belong to a number of clusters that we seek to identify, the performance of the algorithms was measured for a range of tuning parameter values from the set $\left\{2,10,20,50,80,100,200\right\}$. Note that for all algorithms, the parameter value cannot exceed the total number of samples in the data.

For all manifold learning approaches, the following procedure was implemented to tune the optimal parameters of each method per data set:

• 1.

  The method was applied for all parameter values in the set $\left\{2,10,20,50,80,100,200\right\}$.

• 2.

  The K-means algorithm was applied to the low-dimensional embeddings produced for each parameter value

• 3.

  The performance of the chosen method was evaluated quantitatively by computing ACC, NMI, RI and ARI for all tested parameter values.

The optimal parameter value was finally selected based on the evaluation measures. Section “Optimal Parameter Selection” explores how the different approaches are affected by their parameter values. For the other subsections of Section “Result”, the optimal parameter choice per approach was used for the comparison of the multi-view approaches. Section “Optimal Parameter Selection” presents the process of parameter tuning on the synthetic data analysed, and measures the performance of single-view and multi-view manifold learning algorithms. The same process was repeated for the real data analysed (Fig. S4 in Section 4.3 of the Supplemental File).

Synthetic data

A motivational multi-view example was constructed to qualitatively evaluate the performance of multi-view manifold learning algorithms against their corresponding single-view algorithms. Its framework was designed specifically to produce distinct projections of the samples from each data-view. Additional synthetic data sets were generated to explore how the algorithms behave when the separation between the clusters exists, but it is not as explicit as in the motivational example.

All synthetic data were generated using the following process. For the same set of samples, a specified number of data-views were generated, with each data-view capturing different information of the samples. Each data-view, $m$ follows a multivariate normal distribution with mean vector ${\text{μ}}_{\mathbf{m}}=({\mu }_1,\cdots ,{\mu }_{p_m}{)}^T$ and covariance matrix ${\mathrm{\Sigma }}_m=I_{p_m}$, where $p_m$ is the number of features in the $m^{th}$ data-view. The matrix $I_{p_m}$ represents a $p_m\times p_m$ identity matrix. For each data-view, different ${\text{μ}}_{\mathbf{m}}$ values were chosen to distinguish the clusters. Noise, ε, following a multivariate normal distribution with mean ${\mu }_{\epsilon }$ and covariance matrix ${\mathrm{\Sigma }}_{\epsilon }=I_{Pm}$ was added to increase randomness within each data-view. Noise, $\epsilon$, increases the variability within a given data-view. The purpose of this additional variability is to assess whether the algorithms are able to equally capture information from all data-views and are not biased towards the data-view(s) with a higher variability. Thus, noise, $\epsilon$, was only included in selected data-views and opted out from the rest. Although this strategy is equivalent to sampling once using larger variance, the extra noise explicitly distinguishes the highly variable data-view with the rest.

In other words, $X\sim MVN({\mu }_m,{\mathrm{\Sigma }}_m)+\epsilon$, where MV N represents multivariate normal distribution. Distinct polynomial functions (e.g. $h(x)=x^4+3x^2+5$) were randomly generated for each cluster in all data-views and applied to the relevant samples to express non-linearity and to distinguish them into clusters (e.g. the same polynomial function was used for all samples that belong in the same cluster). The last step was performed to ensure that linear dimensionality reduction methods (e.g. PCA) would not successfully cluster the data.

The three synthetic data sets with their characteristics are described next.

Motivational multi-view data scenario

Assume that the truth underlying structure of the data separates the samples into three true clusters as presented in Fig. 1. Each synthetic data-view describes the samples differently, which results in three distinct clusterings, none of which reflects the global underlying truth. In particular, the first view separates only cluster C from the others (View 1 in Fig. 1), the second view separates only cluster B (View 2) and the third view separates only cluster A (View 3). In this scenario, only the third data-view contained an extra noise parameter, $\epsilon$, suggesting a higher variability than the other two data-views.

More clusters than data-views scenario

A synthetic data set that was generated similarly to MMDS but with five true underlying clusters instead of three. The true underlying structure of the each data-view is shown in Fig. 2. In this data set, $p_v=100,\mathrm{\forall }v$ features were generated on $n=500$, equally balanced data samples. In comparison with MMDS, more clusters than data-views scenario (MCS) contains more clusters, but the same number of data-views. Similarly to MMDS and NDS, in this scenario, only the third data-view contained an extra noise parameter, $\epsilon$, suggesting a higher variability than the other two data-views.

Real data

The three real data sets analysed in the study are described below.

Comparison between single-view and multi-view visualisations

Visualising multi-view data can be trivially achieved either by looking at the visualisations produced by each data-view, or by concatenating all features into a long vector. T-SNE, LLE and ISOMAP applied on every single data-view of the MMDS data set separately capture the correct local underlying structure of the respective data-view (Fig. 3). However, by design, they cannot capture the global structure of the data. $\mathrm{S}\mathrm{N}{\mathrm{E}}_{concat}$, $\mathrm{L}\mathrm{L}{\mathrm{E}}_{concat}$ and $\mathrm{I}\mathrm{S}\mathrm{O}\mathrm{M}\mathrm{A}{\mathrm{P}}_{concat}$ represent the trivial solutions of concatenating the features of all data-views before applying t-SNE, LLE and ISOMAP, respectively. These trivial solutions capture mostly the structure of the third data-view, because that data-view has a higher variability between the clusters than the other two.

Multi-SNE, multi-LLE and multi-ISOMAP produced the best visualisations out of all SNE-based, LLE-based and ISOMAP-based approaches, respectively. These solutions were able to separate clearly the three true clusters, with multi-SNE showing the clearest separation between them. Even though m-SNE separates the samples according to their corresponding clusters, this separation would not be recognisable if the true labels were unknown, as the clusters are not sufficiently separated. The visualisation by m-LLE was similar to the ones produced by single-view solutions on concatenated features, while m-ISOMAP bundles all samples into a single cluster.

By visualising the MMDS data set via both single-view and multi-view clustering approaches, multi-SNE has shown the most promising results (Fig. 3). We have shown that single-view analyses may lead to conflicting results, while multi-view approaches are able to capture the true underlying structure of the synthetic MMDS.

Multi-view manifold learning for clustering

It is very common in studies to utilise the visualisation of data to identify any underlying patterns or clusters within the data samples. Here, it is illustrated how the multi-view approaches can be used to identify such clusters. To quantify the visualisation of the data, we applied the K-means algorithm on the low-dimensional embeddings produced by multi-view manifold learning algorithms. If the two-dimensional embeddings can separate the data points to their respective clusters quantitatively with high accuracy via a clustering algorithm, then those clusters are expected to be qualitatively separated and visually shown in two dimensions. For all examined data sets (synthetic and real), the number of clusters (ground truth) within the samples is known, which attracts the implementation of K-means over alternative clustering algorithms. The number of clusters was used as the input parameter, K, of the K-means algorithm and by computing the clustering measures we evaluated whether the correct sample allocations were made.

The proposed multi-SNE, multi-LLE, and multi-ISOMAP approaches were found to outperform their competitive multi-view extensions (m-SNE, m-LLE, m-ISOMAP) as well as their concatenated versions ( $\mathrm{S}\mathrm{N}{\mathrm{E}}_{concat}$, $\mathrm{L}\mathrm{L}{\mathrm{E}}_{concat}$, $\mathrm{I}\mathrm{S}\mathrm{O}\mathrm{M}\mathrm{A}{\mathrm{P}}_{concat}$) (Tables 2 and 3). For the majority of the data sets the multi-SNE approach was found to overall outperform all other approaches.

Figure 4 shows a comparison between the true clusters of the handwritten digits data set and the clusters identified by K-means. The clusters reflecting the digits 6 and 9 are clustered together, but all remaining clusters are well separated and agree with the truth.

Multi-SNE applied on caltech7 produces a good visualisation, with clusters A and B being clearly separated from the rest (Fig. 5B). Clusters C and G are also well-separated, but the remaining three clusters are bundled together. Applying K-means to that low-dimensional embedding does not capture the true structure of the data (Table 2). It provides a solution with all clusters being equally sized (Fig. 5A) and thus its quantitative evaluation is misleading. Motivated by this result, we have further explored the performance of proposed approaches on a balanced version of the caltech7 data set (generated as described in the Section “Real Data”).

Similarly to the visualisation of the original data set, the visualisation of the balanced caltech7 data set shows clusters A, B, C and G to be well-separated, while the remaining are still bundled together (Figs. 5C and 5D).

Through the conducted work, it was shown that the multi-view approaches proposed in the manuscript generate low-dimensional embeddings that can be used as input features in a clustering algorithm (as for example the K-means algorithm) for identifying clusters that exist within the data set. We have illustrated that the proposed approaches outperform existing multi-view approaches and the visualisations produced by multi-SNE are very close to the ground truth of the data sets.

Alternative clustering algorithms, that do not require the number of clusters as input, can be considered as well. For example, Density-based spatial clustering of applications with noise (DBSCAN) measures the density around each data point and does not require the true number of clusters as input (Ester et al., 1996). In situations, where the true number of clusters is unknown, DBSCAN would be preferable over K-means. For completeness of our work, DBSCAN was applied on two of the real data sets explored, with similar results observed as the ones with K-means. The proposed multi-SNE approach was the best-performing method of partitioning the data samples. The analysis using DBSCAN can be found in Figs. S9 and S10 of Section 4.8 in the Supplemental File.

An important observation made was that caution needs to be taken when data sets with imbalanced clusters are analysed as the quantitative performance of the approaches on such data sets is not very robust.

Optimal parameter selection

SNE, LLE and ISOMAP depend on a parameter that requires tuning. Even though the parameter is defined differently in each algorithm, for all three algorithms it is related to the nearest number of global neighbours. As described earlier, the optimal parameter was found by comparing the performance of the methods on a range of parameter values, $S=\left\{2,10,20,50,80,100,200\right\}$. In this section, the synthetic data sets, NDS and MCS, were analysed, because both data sets separate the samples into known clusters by design and evaluation via clustering measures would be appropriate.

To find the optimal parameter value, the performance of the algorithms was evaluated by applying K-means on the low-dimensional embeddings and comparing the resulting clusterings against the truth. Once the optimal parameter was found, we confirmed that the clusters were visually separated by manually looking at the two-dimensional embeddings. Since the data in NDS and MCS are perfectly balanced and were generated for clustering, this approach can effectively evaluate the data visualisations.

On NDS, single-view SNE, LLE and ISOMAP algorithms produced a misclustering error of $0.3$, meaning that a third of the samples was incorrectly clustered (Fig. 6B). This observation shows that single-view methods capture the true local underlying structure of each synthetic data-view. The only exception for NDS is the fourth data-view, for which the error is closer to $0.6$, i.e. randomly assigns the clusters (which follows the simulation design, as it was designed to be a random data-view). After concatenating the features of all data-views, the performance of single-view approaches remains poor (Fig. 6A). The variance of the misclustering error on this solution is much greater, suggesting that single-view manifold learning algorithms on concatenated data are not robust and thus not reliable. Increasing the noise level (either by incorporating additional noisy data-views, or by increasing the dimensions of the noisy data-view) in this synthetic data set had little effect on the overall performance of the multi-view approaches (Table S2, Figs. S5, S6 in Section 4.4. of the Supplemental File).

On both NDS and MCS, multi-LLE and multi-SNE were found to be sensitive to the choice of their corresponding parameter value (Figs. 6A and 7). While multi-LLE performed the best when the number of nearest neighbours was low, multi-SNE provided better results as perplexity was increasing. On the other hand, multi-ISOMAP had the highest NMI value when the parameter was high.

Overall, ISOMAP-based multi-view algorithms were found to be more sensitive to their tuning parameter (number of nearest neighbours). This observation is justified by the higher variance observed in the misclustering error on the optimal parameter value throughout the simulation runs (Fig. 6B). The performance of ISOMAP-based methods improved as the parameter value increased (Fig. 7). However, they were outperformed by multi-LLE and multi-SNE for both synthetic data sets.

Out of the three manifold learning foundations, LLE-based approaches mostly depend on their parameter value to produce the optimal outcome. Specifically, their performance dropped when the parameter value was between $20$ and $100$ (Fig. 6A). When the number of nearest neighbours was set to be greater than $100$ their performance started to improve. Out of all LLE-based algorithms, the highest NMI and lowest misclustering error was obtained by multi-LLE (Figs. 6 and 7). Our observations on the tuning parameters of LLE-based approaches are in agreement with earlier studies (Karbauskaitė, Kurasova & Dzemyda, 2007; Valencia-Aguirre et al., 2009). Both Karbauskaitė, Kurasova & Dzemyda (2007) and Valencia-Aguirre et al. (2009) found that LLE performs best with low nearest number of neighbours and their conclusions reflect the performance of multi-LLE; best performed on low values of the tuning parameter. Even though m-SNE performed better than single-view methods in terms of both clustering and error variability, multi-SNE produced the best results (Figs. 6 and 7). In particular, multi-SNE outperformed all algorithms presented in this study on both NDS and MCS. Even though it performed poorly for low perplexity values, its performance improved for $Perp\ge 20$. Multi-SNE was the algorithm with the lowest error variance, making it a robust and preferable solution.

The four implemented measures (Accuracy, NMI, RI and ARI) use the true clusters of the samples to evaluate the clustering performance. In situations where cluster allocation is unknown, alternative clustering evaluation measures can be used, such as the Silhouette score (Rousseeuw, 1987). The Silhouette score in contrast to the other measures does not require as input the cluster allocation and is a widely used approach for identifying the best number of clusters and clustering allocation in an unsupervised setting.

Evaluating the clustering performance of the methods via the Silhouette score agrees with the other four evaluation measures, with multi-SNE producing the highest value out of all multi-view manifold learning solutions. The Silhouette score of all methods applied on the MCS data set can be found in Fig. S8 of Section 4.7 in the Supplemental File.

The same process of parameter tuning was implemented for the real data sets and their performance is presented in Fig. S4 of Section 4.3 in the Supplemental File. In contrast to the synthetic data, multi-SNE on cancer types data performed the best at low perplexity values. For the remaining data sets, its performance was stable for all parameter values. With the exception of cancer types data, the performance of LLE-based solutions follows their behaviour on synthetic data.

Optimal number of data-views

It is common to think that more information would lead to better results, and in theory that should be the case. However, in practice that is not always true (Kumar, Rai & Hal, 2011). Using the cancer types data set, we explored whether the visualisations and clusterings are improved if all or a subset of the data-views are used. With three available data-views, we implemented a multi-view visualisation on three combinations of two data-views and a single combination of three data-views.

The genomics data-view provides a reasonably good separation of the three cancer types, whereas miRNA data-view fails in this task, as it provides a visualisation that reflects random noise (first column of plots in Fig. 8). This observation is validated quantitatively by evaluating the produced t-SNE embeddings (Table S3 in Section 4.5 of the Supplemental File). Concatenating features from the different data-views before implementing t-SNE does not improve the final outcome of the algorithm, regardless of the data-view combination.

Overall, multi-view manifold learning algorithms have improved the data visualisation to a great extent. When all three data-views are considered, both multi-SNE and m-SNE provide a good separation of the clusters (Fig. 8). However, the true cancer types can be identified perfectly when the miRNA data-view is discarded. In other words, the optimal solution in this data set is obtained when only genomics and epigenomics data-views are used. That is because miRNA data-view contains little information about the cancer types and adds random noise, which makes the task of separating the data points more difficult.

This observation was also noted between the visualisations of MMDS and NDS (Fig. 9). The only difference between the two synthetic data sets is the additional noisy data-view in NDS. Even though NDS separates the samples to their corresponding clusters, the separation is not as clear as it is in the projection of MMDS via multi-SNE. In agreement with the exploration of the cancer types data set, it is favourable to discard any noisy data-views in the implementation of multi-view manifold learning approaches.

It is not always a good idea to include all available data-views in multi-view manifold learning algorithms; some data-views may provide noise which would result in a worse visualisation than discarding those data-views entirely. The noise of a data-view with unknown labels may be viewed in a single-view t-SNE plot (all data-points in a single cluster), or identified, if possible, via quantification measures such as signal-to-noise ratio.

Multi-SNE variations

This section presents two alternative variations of multi-SNE, including automatic weight adjustments and multi-CCA as a pre-training step for reducing the dimensions of the input data-views.

Automated weight adjustments

A simple weight-updating approach is proposed based on the KL-divergence measure from each data-view. This simple weight-updating approach guarantees that more weight is given to the data-views producing lower KL-divergence measures and that no data-view is being completely discarded from the algorithm.

Recall that $KL(P||Q)\in [0,\mathrm{\infty })$, with $KL(P||Q)=0$, if the two distributions, P and Q, are perfectly matched. Let $\mathbf{k}=(k^{(1)},\cdots ,k^{(M)})$ be a vector, where $k^{(m)}=KL(P^{(m)}||Q),\mathrm{\forall }m=\{1,\cdots ,M\}$ and initialise the weight vector $\mathbf{w}=(w^{(1)},\cdots ,w^{(M)})$ by $w^{(m)}=\frac{1}{M},\mathrm{\forall }m$. To adjust the weights of each data-view, the following steps are performed at each iteration:

• 1.

  Normalise KL-divergence by $k^{(m)}=\frac{k^{(m)}}{{\sum }_i^M{k}^{(i)}}$. This step ensures that $k^{(m)}\in [0,1],\mathrm{\forall }m$ and that ${\sum }_m{k}^{(m)}=1$.

• 2.

  Measure the weights for each data-view by $w^{(m)}=1-k^{(m)}$. This step ensures that the data-view with the lowest KL-divergence value receives the highest weight.

Based on the analysis in Section “Optimal Number of Data-views”, we know that cancer types and NDS data sets contain noisy data-views and thus multi-SNE performs better when they are entirely discarded. Here, we assume that this information is unknown and the proposed weight-updating approach is implemented on those two data sets to test if the weights are being adjusted correctly according to the noise level of each data-view.

The proposed weight-adjustment process, which looks at the produced KL-divergence between each data-view and the low-dimensional embeddings, distinguishes which data-views contain the most noise and the weight values are updated accordingly (Fig. 10B). In cancer types, transcriptomics (miRNA) receives the lowest weight, while genomics (Genes) was given the highest value. This weight adjustment comes in agreement with the qualitative (t-SNE plots) and quantitative (clustering) evaluations performed in Section “Optimal Number of Data-views”. In NDS, $X^{(4)}$ which represents the noisy data-view received the lowest weight, and the other data-views had around the same weight value, as they all impact the final outcome equally.

The proposed weight-adjustment process updates the weights at each iteration. For the first $100$ iterations, the weights are not changing, as the algorithm adjusts to the produced low-dimensional embeddings (Fig. 10B). In NDS, the weights converge after $250$ iterations, while in cancer types, they are still being updated even after $1,000$ iterations. The changes recorded are small and the weights can be said to have stabilised.

The low-dimensional embeddings produced in NDS with weight adjustments separate clearly the three clusters, an observation missed without the implementation of the weight-updating approach (Fig. 10A); it resembles the MMDS (i.e. without noisy data-view) multi-SNE plot (Fig. 9). The automatic weight-adjustment process identifies the informative data-views, by allocating them a higher weight value than to the noisy data-views. This observation was found to be true even when a dataset contains more noise than informative data-views (Table S2, Figs. S5, S6 in Section 4.4 of the Supplemental File).

The embeddings produced using multi-SNE with the automated weight adjustments in cancer types do not separate the three clusters as clearly as multi-SNE without the noisy data-view. A more clear separation is obtained compared to multi-SNE on the complete data set without weight adjustments.

The weights produced by this weight adjustment approach can indicate the importance of each data-view in the final lower-dimensional embedding. For example, data-views with very low weights may be assumed futile and a better visualisation may be produced if those data-views are discarded. The actual weights assigned to each data-view do not have any further meaning.

Multi-CCA as pre-training

As mentioned earlier, van der Maaten & Hinton (2008) proposed the implementation of PCA as a pre-training step for t-SNE to reduce the computational costs, provided that the fraction of variance explained by the principal components is high. In this article, pre-training via PCA was implemented in all variations of SNE. Alternative linear dimensionality reduction methods may be considered, especially for multi-view data. In addition to reducing the dimensions of the original data, such methods can capture information between the data-views. For example, canonical correlation analysis (CCA) captures relationships between the features of two data-views by producing two latent low-dimensional embeddings (canonical vectors) that are maximally correlated between them (Hotelling, 1936; Rodosthenous, Shahrezaei & Evangelou, 2020). Rodosthenous, Shahrezaei & Evangelou (2020) demonstrated that multi-CCA, an extension of CCA that analyses multiple (more than two) data-views, would be preferable as it reduces over-fitting.

This section demonstrates the application of multi-CCA as pre-training in replacement of PCA. This alteration of the multi-SNE algorithm was implemented on the handwritten digits data set. Multi-CCA was applied on all data-views, with six canonical vectors produced for each data-view (in this particular data set $min(p_1,p_2,p_3,p_4,p_5,p_6)=6$). The variation of multi-CCA proposed by Witten & Tibshirani (2009) was used for the production of the canonical vectors, as it is computationally cheaper compared to others (Rodosthenous, Shahrezaei & Evangelou, 2020). By using these vectors as input features, multi-SNE produced a qualitatively better visualisation than using the principal components as input features (Fig. 11). By using an integrative algorithm as pre-training, all $10$ clusters are clearly separated, including six and nine. Quantitatively, clustering via K-Means was evaluated with ACC = 0.914, NMI = 0.838, RI = 0.968, ARI = 0.824. This evaluation suggests that quantitatively, it performed better than the $10$-dimensional embeddings produced multi-SNE with PCA as pre-training.