We recall some geometric properties of simplices, which our method is based on. For a set of n points (x ₁, … , x _n ), the associated n-simplex is the polytope of vertices (x ₁, … , x _n ) (a 3-simplex is a triangle, a 4-simplex is a tetrahedron and so on). The height h(V _n , x) of a point x belonging to a n-simplex V _n can be obtained as (Sommerville, 1929) h V n , x = n V n V n − 1 , (1) where V _n is the volume of the n-simplex, and V _n−1 is the volume of the (n − 1)-simplex obtained by removing the point x. V _n and V _n−1 can be computed using the pairwise distances only, with the Cayley-Menger formula (Sommerville, 1929): V n = | det C M n | 2 n ⋅ n ! 2 , (2) where det(CM _n ) is the determinant of the Cayley-Menger matrix CM _n , that contains the pairwise distances d i , j = x i − x j , as C M n = 0 1 1 … 1 1 1 0 d 1,2 2 … d 1 , n 2 d 1 , n + 1 2 1 d 2,1 2 0 … d 2 , n 2 d 2 , n + 1 2 … … … … … … 1 d n , 1 2 d n , 2 2 … 0 d n , n + 1 2 1 d n + 1,1 2 d n + 1,2 2 … d n + 1 , n 2 0 . (3)

We now consider a dataset X of size N × d, where N is the sample size and d the dimension of the data. We associate with X a matrix D of size N × N, which represents all the pairwise distances between observations of X. We also assume that the data points can be mapped into a vector space with regular observations that form a main subspace of unknown dimension d* with some small noise, and additional orthogonal outliers of relatively large orthogonal distance to the main subspace (Figure 1A). Our proposed method aims to infer from D the dimension of the main data subspace d*, using the geometric properties of simplices with respect to their number of vertices: Consider a (n + 2)-simplex containing a data point x _i and its associated height, that can be computed using Eq. 1 in Section 2.1. When n < d* and for S large enough, the distribution of heights obtained from different simplices containing x _i remains similar, whether x _i is an orthogonal outlier or a regular observation (see Figure 1B). In contrast, when n ≥ d*, the median of these heights approximately yields the distance of x _i to the main subspace (Figure 1C). This distance should be significantly larger when x _i is an orthogonal outlier, compared with regular points, for which these distances are tantamount to the noise.

To estimate d* and for a given dimension n tested, we thus randomly sample, for every x _i in X, S(n + 2)-simplices containing x _i , and compute the median of the heights h i n associated with these S simplices. Upon considering, as a function of the dimension n tested, the distribution of median heights ( h 1 n , … , h N n ) (with 1 ≤ i ≤ N), we then identify d* as the dimension at which this function presents a sharp transition towards a highly peaked distribution at zero. To do so, we compute h ̃ n , as the mean of ( h 1 n , … , h N n ) , and estimate d* as n ̄ = argmax n h ̃ n − 1 h ̃ n . (4)

Furthermore, we detect orthogonal outliers using the distribution obtained in n ̄ , as the points for which h i n ̄ largely stands out from h ̃ n ̄ . To do so, we compute σ n ̄ the standard deviation observed for the distribution ( h 1 n ̄ , … , h N n ̄ ) , and obtain the set of orthogonal outliers O as O = i | h i n ̄ > h ̃ n ̄ + c × σ n ̄ , (5) where c > 0 is a parameter set to achieve a reasonable trade-off between outlier detection and false detection of noisy observations. Our implementation uses c = 3 by default (following the three σ rule Pukelsheim (1994), and which corresponds to ∼ 99.9 % of a Gaussian distribution being conserved), value which was also used in our experiments. In case users possess prior information or want to control the fraction of detected outliers, the value of c may be modified, with increasing c making the detection stricter. Also note that our method introduces another parameter S, as it samples S simplices to calculate the median of the corresponding heights. Therefore, S should be large enough so the resulting sample median well approximates the global median. Assuming the heights being sampled from a continuous distribution, this can be guaranteed as the sample median is asymptotically normal, with mean equal to the true median and the standard deviation proportional to 1 S (Rider, 1960).

The method presented in the previous section assumes that at dimension d*, the median height calculated for each point reflects the distance to the main subspace. This assumption is valid when the fraction of orthogonal outliers is small enough, so that the sampled n-simplex likely contains regular observations only, aside from the evaluated point. However, if the number of outliers gets large enough so that a significant fraction of n-simplices also contains outliers, then the calculated heights would yield the distance between x _i and an outlier-containing hyperplane, whose dimension is larger than a hyperplane containing only regular observations. The apparent dimensionality of the main subspace would thus increase and generates a positive bias on the estimate of d*.

Specifically, if X contains a fraction of p outliers, and if we consider o _n,p,N the number of outliers drawn after uniformly sampling n + 1 points (to test the dimension n), then o _n,p,N follows a hypergeometric law, with parameters n + 1, the fraction of outliers p = N _o /N, and N. Thus, the expected number of outliers drawn from a sampled simplex is (n + 1) × p. After estimating n ̄ (from Section 2.2), and finding a proportion of outliers p ̄ = | O | / N using Eq. 5, we hence correct n ̄ a posteriori by substracting the estimated bias δ, as the integer part of the expectation of o _n,p,N , so the debiased dimensionality estimate n* is n * = n ̄ − ⌊ n ̄ + 1 × p ⌋ . (6)

Upon identifying the main subspace containing regular points, our procedure finally corrects the pairwise distances that contain outliers in the matrix D, in order to apply a MDS that projects the outliers in the main subspace. In the case where the original coordinates cannot be used (e.g., as a result of some transformation or if the distance is non Euclidean), we perform the two following steps: 1) We first apply a MDS on D to place the points in a euclidean space of dimension d, as a new matrix of coordinates X ̃ . 2) We run a PCA on the full coordinates of the estimated set of regular data points (i.e., X ̃ \ O ), and project the outliers along the first n ̄ * principal components of the PCA, since these components are sufficient to generate the main subspace. Using the projected outliers, we accordingly update the pairwise distances in D to obtain the corrected distance matrix D*. Note that in the case where D derives from a euclidean distance between the original coordinates, we can skip step 1), and directly run step 2) on the full coordinates of the estimated set of regular data points.

The overall procedure, called DeCOr-MDS, is described in Algorithm 1. The values for the parameters S and c were set by default and in our experiments to S = 100 and c = 3. We also provide an implementation in Python 3.8.10 available on this github repository: https://github.com/wxli0/DeCOr-MDS.

Algorithm 1

DeCOr-MDS.

Input D the pairwise distance matrix of the dataset of size N × d, E _dim the set of dimensions ( ≤ d ) to be tested, c and S user-specified constants

The complexity of the algorithm can be briefly evaluated as follows.

1. Given one n-simplex, the volume computation has a complexity of O ( n 3 ) . Since we compute the height for S simplices and repeat the process for all E _dim dimensions, the total complexity of this step amounts to O ( S N n 3 E dim ) ,

Note that with the tested dimensions being smaller than the data dimension (n, E _dim < d), and the number of simplices being significantly smaller than the total number of data points (S ≪ N), the burden of evaluating the simplices (step 1) and correcting outliers (step 4) is in practice less than the cost of the PCA (step 2) and MDS (step 3).

We propose a robust method to reduce and infer the dimensionality of a dataset from its pairwise distance matrix, by detecting and correcting orthogonal outliers. The method, called DeCOr-MDS, can be divided into three sub-procedures detailed in Sections 2.2–2.4, with the overall algorithm provided in Section 2.5. The first procedure detects orthogonal outliers and estimates the subspace dimension using the statistics of simplices that are sampled from the data, using Eqs 4, 5. The second procedure corrects for potential bias in estimated dimension when the fraction of outlier is large. The third procedure corrects the pairwise distance of the original data, by replacing the distance to orthogonal outliers by that to their estimated projection on the main subspace. In the next sections, we report the results obtained upon running the procedure on synthetic and various biological datasets, that demonstrate the performance and accuracy of the method. For all these experiments, we also reported the runtime in Supplementary Table S1, showing how the method can be used in practice with reasonable time on experimental datasets (less than 10 min in our workstation, with x86_64 CPU, 132 GB RAM and 447 GB disk storage).

We first illustrate and evaluate the performance of the method on synthetic datasets, (for a detailed description of the datasets and their generation, see the Methods Section 2.6). On a simple dataset of points forming a 2D cross embedded in 3D (Figure 2A), we observed that the MDS is sensitive to the presence of orthogonal outliers and distorts the cross when reducing the data in 2D (Figure 2B). In contrast, our procedure recovers the original geometry of the uncontaminated dataset, with the outliers being correctly projected (Figure 2C). The same results were obtained when sampling regular points from a 2D plane (Supplementary Figure S1). We further tested higher dimensions, and illustrate in Figure 3A how the distribution of heights becomes concentrated around 0, when testing for the true dimension (d* = 10), as suggested in the Methods Section 2.2. As a result, our method allows to infer the main subspace dimension from Eq. 4, as shown in Figure 3B. In addition, the procedure accurately corrects the pairwise distances to orthogonal points with the distances to their projections on the main subspace, as shown in Figure 3C.

When the dimension of the subspace and fraction of outliers get significantly large, we also illustrate the importance of the correction step (see Methods Section 2.3), due to the sampling of simplices that contain several outliers. Upon using synthetic datasets with d* = 40 and varying the number (fraction) of outliers from 20 (2%) to 100 (10%), we observe this bias appearing before correction, with d* being overestimated by 2 or 3 dimensions (Figure 4). Using the debiased estimate n* from Eq. 6 successfully reduced the bias, with an error ≤ 1 for all the parameters tested.

We further show how DeCOr-MDS can be broadly applied to biological data, ranging from images to high throughput sequencing. We first studied a dataset of single cell images, from osteosarcoma cells (see Figure 5A), which were processed to extract from their contour a 100 × 2 array of xy coordinates representing a discretization of a closed curve (see Dataset Section 2.6). We obtained a pairwise distance matrix on this set of curves by using the so-called Square Root Velocity (SRV) metric, which defines a Euclidean distance on the space of velocities that derive from a regular parameterization of the curve (Srivastava et al., 2010; Miolane et al., 2020). Using DeCOr-MDS, we found a main subspace of dimension 2 (Figure 5B), with 14 (7%) outliers detected among the 207 cells of this dataset. The comparison between the resulting embedding and that obtained from a simple MDS is shown in Supplementary Figure S2, and reveals that outliers, when uncorrected, affect the embedding coordinates, while our correction mitigates it. By examining in more details the regular and inferred outlier cells (Figure 5C, with all cell shapes shown in Supplementary Figure S3), we found regular observations to approximately describe elliptic shapes, which is in agreement with the dimension found, since ellipses are defined by 2 parameters. One can also visually interpret the orthogonal outliers detected as being more irregular, with the presence of more spikes and small protusions. Interestingly, the procedure also identified as outliers some images containing errors, due to bad cropping or segmentation (with 2 cells shown instead of one), which should thus be removed of the dataset for downstream analysis.

As another example of application to biological data, we next considered a dataset from the Human Microbiome Project (HMP). The Human Microbiome Project aims at describing and studying the microbial contribution to the human body. In particular, genes contributed by microbes in the gut are of primary importance in health and disease (Turnbaugh et al., 2007). The resulting data is an array which typically contains the abundance of different elements of the microbiome (typically 10²–10³), denoted phylotypes, measured in different human subjects. To analyze such high dimensional datasets, dimensionality reduction methods including MDS (often denoted Principal Coordinates Analysis PCoA), are typically applied and used to visualize the data (Brooks et al., 2018; Trevelline and Kohl, 2022; Zhou et al., 2022).

To assess our method incrementally, we restricted first the analysis to a representative specific site (nose), yielding a 136 × 425 array that was further normalized to generate Euclidean pairwise distance matrices (see Material and Methods Section 2.6 for more details). Upon running DeCOr-MDS, we estimated the main dimension to be 3, with 9 (6.62%) orthogonal outliers detected, as shown in Figure 6A. This is also supported by another study that the estimated dimension of HMP dataset is 2 or 3 (Tomassi et al., 2021). We also computed the average distance between these orthogonal outliers and the barycenter of regular points in the reduced subspace, and obtained a decrease from 1.21 when using MDS to 0.91 when using DeCOr-MDS. This decrease suggests that orthogonal outliers get corrected and projected closer to the regular points, to improve the visualization of the data in the reduced subspace, like in our experiments with the synthetic datasets (Figure 2 and Supplementary Figure S2). In Figure 6B, we next aggregated data points from another site (throat) to study how the method performs in this case, yielding a 270 × 425 array that was further normalized to generate Euclidean pairwise distance matrices. As augmenting the dataset brings a separate cluster of data points, the dimension of the main dataset was then estimated to be 2, with 13 (5%) orthogonal outliers detected, as shown in Figure 6B. The average distance between the projected outliers and the barycenter of projected regular points are approximately the same when using MDS (1.46) as when using DeCOr-MDS (1.45) for nose, and are also approximately the same when using MDS (1.75) to when using DeCOr-MDS (1.74) for throat. This decrease also suggests that orthogonal outliers get corrected and projected closer to the regular points.

We further evaluated DeCOr-MDS on single cell RNA-seq (scRNA-seq) data. In general, analyzing scRNA seq data requires dimensionality reduction for visualization (including MDS-based methods Canzar et al. (2021); Senabouth et al. (2019)), and specific quality control procedures to mitigate various technical artifacts McCarthy et al. (2017); Luecken and Theis (2019). We applied our method as a potentially relevant tool for this purpose. We applied DeCOr-MDS first on a dataset containing the expression level of 1063 cells for 2,601 genes (detailed in Methods Section 2.6). We found the dimension of the main subspace to be 3, with 77 (7%) outliers detected. In Figures 7A,B, we compared the embeddings in 3D using MDS and DeCOr-MDS. Similarly to the previous experiments, the mean distance between the orthogonal outliers and barycenter of regular points in the reduced subspace decreases when using DeCOr-MDS (from 4.22 to 2.51), improving the visualization of regular points. In Figure 7C, we further examined the drop-out rates (indicating zero count for a given gene) of the cells among the top 500highly expressed genes, determined by the median of counts per gene. Among these highly expressed genes, we identified 97.4% of the detected outliers that have drop-out rates greater than 0.95, while this was the case for 27.4% of the regular cells. Upon performing a pairwise t-test on the total counts for the top 500 highly expressed genes from the outlier group and the regular cell group, we found that the total counts are significantly different between the two groups (p-value < 0.001). Therefore, our method led to detect some outliers associated with high drop-out counts for highly expressed genes, which were not captured by the standard processing and quality control methods used in the first place.