Sailing in high-dimensional spaces: Low-dimensional embeddings through angle preservation

## Abstract

Abstract Low-dimensional embeddings (LDEs) of high-dimensional data are ubiquitous in science and engineering. They allow us to quickly understand the main properties of the data, identify outliers and processing errors, and inform the next steps of data analysis.
As such, LDEs have to be faithful to the original high-dimensional data, i.e., they should represent the relationships that are encoded in the data, both at a local as well as global scale.
The current generation of LDE approaches focus on reconstructing local distances between any pair of samples correctly, often outperforming traditional approaches aiming at all distances.
For these approaches, global relationships are, however, usually strongly distorted, often argued to be an inherent trade-off between local and global structure learning for embeddings. We suggest a new perspective on LDE learning, reconstructing angles between data points.
We show that this approach, Mercat , yields good reconstruction across a diverse set of experiments and metrics, and preserve structures well across all scales.
Compared to existing work, our approach also has a simple formulation , facilitating future theoretical analysis and algorithmic improvements.

## 1 Introduction

Figure: Figure 1: Visual abstract. Existing work (top) optimizes low-dimensional embeddings to reconstruct distances, focusing on reconstruction of local structures (smaller distances), leading to distortion or breaking of global structures (larger distances). We suggest (bottom) to reconstruct angles between any three points, embedding on the sphere 2D sphere $\mathbb{S}^{2}$, capturing structures at any scale.
Refer to caption: extracted/5666939/figs/mercat_idea.png

A key aspect of modern data analysis is data visualization, which allows domain experts to directly interact with their data and get a glimpse into its hidden patterns. Usually employed early in an analysis, such a visualization can help to identify errors in data recording and preprocessing, discover outliers, and overall structure of the data, informing next processing steps or choice of further analysis.
Typically, however, data is high-dimensional and hence does not lend itself for direct visualization. Methods computing low-dimensional embeddings (LDEs) take on this problem, computing an embedding of the data in typically 2 or 3 dimensions that can be visually perceived by humans. Such methods are nowadays part of the standard data analysis pipelines in e.g. biology where LDEs are used to visualize complex gene regulatory data [], or as an exploratory tool in machine learning, where it is used to investigate latent spaces of neural networks [, , ].

To get a proper understanding of the data and then make informed decisions, the LDEs have to be faithful to the original data: local structures should be perceivable, but also global relationships of these structures should be appropriately reflected.
While several widely used methods to obtain LDEs exist, it has been observed that often they only reconstruct local structures faithfully, while neglecting global structures [, , , , ].
This leads to a loss of information in the embedding space, where for data consisting of clusters most inter-cluster information is lost [], and for data with manifold structures, the manifold gets absurdly distorted or torn, capturing only locally faithful information that make it hard to reason about the data as a whole [, , ] (cf. Fig. 2).

This drawback likely comes by design, as current state-of-the-art approaches mainly focus on the correct reconstruction of local distances, neglecting long-range distances in the process. A common argument for this approach is that it is impossible to compress all information in the high-dimensional data into the low-dimensional space. A theoretical barrier is expected to exist which reflects the fundamental trade-off between local and global structure preservation inherent in LDEs. Due to the iterative nature of the optimization and the rather complex objective functions involved in most of the existing approaches, theoretical understanding along this direction is still very limited [, , ]. Nevertheless, empirical studies [, ] suggest it is unlikely that existing methods have already reached such a critical point, that any further improvement on the global structure preservation has to come with a compromise in the faithfulness of local structure preservation.

In this work, we undertake a new approach that is distinct from the common paradigm of reconstructing (local) distances and propose a simple new LDE approach preserving both local and global structures well.
Inspired by a breakthrough in navigation in the 16th century, the Mercator Projection, we aim for a representation that focuses on the reconstruction of angles.
The motivation for the Mercator Projection—a projection of planet Earth on a 2D map—was that flat lines on the map represent routes of constant bearing, greatly simplifying navigation with a compass. This 2D map, which most of us recognize from world maps in an atlas or online mapping services, implements the idea of preserving angles locally at every point, thus keeping relative orientation of objects (landmasses) intact.
Mathematically, such a map is called a conformal map [, ], where every angle between two curves that cross each other is preserved.
We follow this idea and introduce the concept of a angle-approximating embeddings (Fig. 1,) which we define as a lower dimensional representation of the original data in which angles between data points are preserved, mapping a high-dimensional dataset on the 2D unit sphere $\mathbb{S}^{2}$, which can be directly visualized.

Our approach, Mercat, is both theoretically appealing due to its simplicity, and practically useful.
We show on challenging toy examples, synthetic data studies, and real data sets, that the embeddings are not only visually more faithful to the original data but also quantitatively great in reconstruction of both local and global structure.
It, hence, serves as a basis for new developments in LDE theory and practice.
Concretely, our contributions are (i) we propose a new paradigm for computing low-dimensional embedding by optimizing for reconstruction of angles rather than distances, (ii) provide efficient algorithmic ideas to compute such an LDE in practice, (iii) give empirical and theoretical justifications for our algorithmic ideas, and (iv) provide extensive evaluation on synthetic and real-world data against state-of-the-art approaches with a diverse set of metrics.

## 2 Related Work

We will here focus on reviewing the research on LDEs most relevant to ours, in particular different types of embedding strategies. Perhaps the most widely known dimensionality reduction method is principal component analysis (PCA) [], followed by seminal work on multidimensional scaling [], self-organizing maps [], and Laplacian eigenmaps []. What all of these methods have in common is that their objective usually focuses on getting the larger distances right. With empirical evidence and the methodological insight that data often lies on an intrinsically low-dimensional manifold, following work such as local linear embedding [], Isomap [], and Hessian eigenmaps [] focus on getting local distances right and modeling relationships non-linearly. However, these methods rely on stringent manifold assumptions, restricting their scalability for large and high-dimensional data, and making their practical performance susceptible to noise and data outliers. More recently, a family of low-dimensional embedding algorithms based on ideas of stochastic neighbor embeddings (SNE) [], with t-distributed SNE (tSNE) [] and the closely related Uniform Manifold Approximation (UMAP) [] being its most prominent representatives, have become extremely popular in data analysis and scientific research, especially in the field of molecular biology [, ]. These algorithms again focus on the reconstruction of local neighborhoods, but have been found more scalable and more robust to high-dimensional noisy data sets, compared with previous methods.
Recent research continued to improve this line of research by—for example—improving the computational efficiency using fast-fourier-based optimization techniques for tSNE [], while keeping the original paradigm of tSNE and UMAP largely untouched.
NCVis [] instead formulate a noise-contrastive learning problem, which is fast in practice, but still emphasizes local distances.

While widely employed, both tSNE and UMAP as well as related approaches suffer from several known issues, one of them being that densities are not properly preserved in the embeddings – two vastly differently sized clusters are mapped to the same amount of space in the embedding. Recent works focused on solving these issues [, ].
Other debates focus on the severe distortions of global structures caused by the neglect of long-range distances in the reconstruction [, , ], which however has not reached a convincing solution to this problem, so far leaving open whether this issue is due to a theoretical trade-off inherent in the problem or a problem of modeling choice.

## 3 Angle-preserving low-dimensional embeddings

The key idea of our method is motivated by the central issue of state-of-the-art LDE approaches, which is the poor reconstruction of long-range, or “global”, relationships, such as the orientation of sub-manifolds or the relative locations of clusters (see Fig. 2).
This problem comes by design, as concurrent work, e.g., tSNE and UMAP, focuses on reconstructing local distances well. This approach so far outperforms traditional methods aiming to reconstruct all distances, such as MDS [], revealing more meaningful structure.
Often discussed as a necessary theoretical trade-off between good reconstruction of local versus global relationships, we, however, instead hypothesize that the problem is inherent in the modeling of distances.
Inspired by the Mercator projection, which revolutionized navigation in the 16th century by providing an angle-preserving (conformal) 2D map of the earth, we suggest to compute an embedding that approximately reconstructs angles between any three data points (angles within each triangle, see Fig. 1). As such, it inherently balances global and local relationships by being independent of the scale of the structures.

Following the idea of a map similar to a mercator projection, which is a map of fixed size and is without borders (i.e., “leaving” the map on the left means “entering” it on the right), we compute an embedding on the unit sphere. The unit sphere is a 2-dimensional space of fixed area and without borders, which allows to embed complex and possibly periodic patterns commonly arising from biological applications such as cell linage [], cell cycle [, ] and circadian rhythm [], but at the same time lends itself for efficient computation of angles and (geodesic) distances.

### 3.1 Formal description

For data $X=\{X_{i}\}_{1\leq i\leq n}\subset\mathbb{R}^{d}$ of $n$ samples and $d$ features, we are interested in a good reconstruction of $X$ in a low-dimensional space, particularly the (unit) 2-sphere $\mathbb{S}^{2}$.
As introduced above, we consider an approximation of a conformal embedding for our LDE, i.e. an embedding that approximately reconstructs angles from the data in the high-dimensional space, thus orienting both local as well as global structures properly to each other.
More formally, we are searching for a map $X_{i}\mapsto Y_{i}$, where $Y_{i}\subset\mathbb{S}^{2}$, such that for any sample $i$ and pair of samples $j,k$,
the (Euclidean) angle between $X_{j}$ and $X_{k}$ measured at $X_{i}$ should be reconstructed in $Y=\{Y_{i}\}_{1\leq i\leq n}$, i.e., $\angle X_{j}X_{i}X_{k}\approx\angle Y_{j}Y_{i}Y_{k}$. Here $\angle Y_{j}Y_{i}Y_{k}$ is the angle between the shortest paths $\overline{Y_{i}Y_{j}}$ and $\overline{Y_{i}Y_{k}}$ on the $\mathbb{S}^{2}$ sphere, i.e., the angle between the corresponding geodesics.
Using the definition of Euclidean inner product in $\mathbb{R}^{d}$, we get
$\arccos\left(\frac{(X_{j}-X_{i})^{\top}(X_{k}-X_{i})}{||X_{j}-X_{i}||~{}||X_{k
}-X_{i}||}\right)\approx\angle Y_{j}Y_{i}Y_{k}.$
We parameterize each point $Y_{i}$ in $\mathbb{S}^{2}$ by two parameters $(\phi_{i},\theta_{i})$, which correspond to longitude and latitude on the 2-sphere.
The angles can then be computed in terms of these coordinates using the following relations.

*Side-to-angle formula (spherical law of cosines).*
Let $\Delta ABC$ be a triangle on the sphere, with $a=\overline{BC}$, $b=\overline{AC}$, $c=\overline{AB}$, $\alpha=\angle(CAB)$, $\beta=\angle(CBA)$, and $\gamma=\angle(ACB)$. Then we have
$\cos\alpha=\frac{\cos a-\cos b\cos c}{\sin b\sin c},$ $\cos\beta=\frac{\cos b-\cos c\cos a}{\sin c\sin a},$ $\cos\gamma=\frac{\cos c-\cos a\cos b}{\sin a\sin b}.$

*Vertex-to-side formula.*
For $Y_{i}=(\phi_{i},\theta_{i})$ and $Y_{j}=(\phi_{j},\theta_{j})$ on the sphere, it follows that the geodesic distance between $Y_{i}$ and $Y_{j}$ is
$\overline{Y_{i}Y_{j}}=d(Y_{i},Y_{j})=\sqrt{2-2[\sin\phi_{i}\sin\phi_{j}\cos(
\theta_{i}-\theta_{j})+\cos\phi_{i}\cos\phi_{j}]}.$

Objective. Considering the root mean square deviation of angles in $Y$ from those corresponding in $X$, we get a differentiable objective

$$ $\mathcal{L}(X,Y)=\bigg{(}\frac{1}{n}\sum_{i}\frac{1}{(n^{2}-n)/2}\sum_{{(i,j,k ):k>j\text{ and }j,k\neq i}}||\angle X_{j}X_{i}X_{k}-\angle Y_{j}Y_{i}Y_{k}||^ {2}_{2}\bigg{)}^{1/2},$ (1) $$

which we can optimize as
$\operatorname*{arg\,min}_{Y\subset\mathbb{S}^{2}}\mathcal{L}(X,Y).$
Through the parameterization by longitude and latitude on the 2-sphere as defined above,
we can optimize directly on the sphere by standard gradient descent on $\frac{\partial\mathcal{L}(X,Y)}{\partial\phi},\frac{\partial\mathcal{L}(X,Y)}{
\partial\theta}$.
Having all components of our approach together, we give the pseudocode of our method—named as Mercat in reminiscence of the inspirational idea from the Mercator projection—in Algorithm 1
.

Figure: Algorithm 1 Mercat

### 3.2 Computational and statistical strategies

Initialization.
Low-dimensional embedding techniques greatly benefit from a good initialization []. For a good initial embedding we follow the established strategy and consider the first two principal components as initialization, wrapping them around the sphere.
In particular, let $PC1$ and $PC2$ be the points in $X$ projected onto the first and second principal component, respectively.
We then compute the initial longitudes as $0.6\pi\frac{PC1-\min(PC1)}{\max(PC1)-\min(PC1)}+0.2\pi$ and initial latitudes as $0.6\pi\frac{PC2-\min(PC2)}{\max(PC2)-\min(PC2)}+0.2\pi$.
We thus roughly distribute the points on a half-sphere while keeping the initial estimates away from the poles(^1^11Having points close to a pole leads to slow optimization as the loss landscape is flat around them.).

Two simplifying computational tricks. In practice, we can employ two computational tricks to accelerate the optimization and improve numerical stability.
First, we drop the arc-cosine from angle computations, i.e., we compute differences between normalized dot products, which is a strictly monotone transformation of the original formulation that does not change the optima.
Second, we compute angles on $Y$ by pure linear algebra, using the dot product between normals of two planes (see App. Sec. A.1), which can be much faster when employed on modern hardware such as GPUs.
We further incorporate two statistical techniques, which speed up the algorithm and improve its scalability while making the embeddings robust to noise and high dimensionality.

Angle evaluation after spectral denoising. The first statistical strategy is to denoise the original high-dimensional data using spectral methods before evaluating their Euclidean angles. It is known that in high dimensions, the Euclidean angles between data points can be sensitive to noise perturbations and may suffer severely from the effect of high dimensionality. We argue that in many applications the observed high-dimensional data points are only noisy versions of some latent noiseless samples incorporating certain low-dimensional signal structures. As such, the quantity of interest should be the Euclidean angles among the noiseless samples, of which the angles among the original noisy high-dimensional data can be very poor estimates (see [, , ] and Theorem 2 below). To overcome such limitations, instead of directly calculating the angles among the original high-dimensional data points, we propose to first apply a principal component analysis (PCA) to the data matrix $X$, to obtain denoised low-dimensional spectral embeddings given by the leading $r$ principal components. After that, we use the Euclidean angles calculated from such low-dimensional spectral embeddings to estimate the angles among the noiseless samples. Note that while concurrent work on low-dimensional embeddings often considers a similar approach of projecting high-dimensional data to a few principal components before embedding them, the usual justification of this procedure was based on empirical observations of more stable results. We here provide rigorous theoretical justification of our denoising procedure in Section 3.3.

Subsampling.
In practice, computing all angles in every iteration would incur a computational cost in $O(kn^{3})$ for $k$ iterations and $n$ datapoints.
While much of it can be efficiently computed by using linear algebra instead of trigonometry, it is still hard to scale to large datasets.
We thus investigate whether it is indeed necessary to compute all angles in every iteration.
For an empirical study, we consider a real dataset about single-cell gene expression of human hematopoiesis [], a typical application for low-dimensional embeddings. We sample $n=500$ points and compute for each point $X_{i}$ all angles at that point, i.e., all $\angle X_{j}X_{i}X_{k},j\neq k\neq i$, yielding matrices of cosine-angles $\Theta_{i}[j,k]=\cos(\angle X_{j}X_{i}X_{k})$.
We then compute a singular value decomposition for each of these matrices (cf App. Fig.4a), which show that the matrices have only few large singular values. Furthermore, computing the effective rank of the matrix [], which give an estimate of the intrinsic dimensionality of the matrix based on the singular values, we observe that the effective rank is very low (cf App. Fig.4b), with a mean of $13.8$.
Based on these insights, rather than computing all angles at every point, we suggest to sample a fraction of points at random each time we compute angles at point $X_{i}$, i.e., we consider $\angle X_{j}X_{i}X_{k},j\neq k$ and $j,k\in S(n)\setminus\{i\}$, where $S(n)$ is a random subset of $[n]$.
For the remainder of the paper, in every iteration for each point $X_{i}$ we will draw 64 other points uniformly at random and compute angles at $\angle X_{j_{1}}X_{i}X_{j_{2}}$, where $j_{1},j_{2}$ are from these sampled subsets,
effectively reducing the computational costs to $\mathcal{O}(kn)$.

### 3.3 Theoretical justifications

Here we provide theoretical justification for the spectral estimation of the true Euclidean angles of the input data. For a theoretical reasoning of the efficacy of angle subsampling, we refer to App. Sec. B.3.
A given input data is usually noisy and high-dimensional but contains low-dimensional structures. To fix ideas, we first introduce the statistical framework of the spiked population model [, , , , , ]. We assume the high-dimensional data
matrix $X\in\mathbb{R}^{n\times d}$ satisfies
$X=\sum_{i=1}^{r}\sqrt{\lambda_{i}}\mathbb{u}_{i}\mathbf{y}_{i}^{\top}+Z=UY^{
\top}+Z$
where $Z\in\mathbb{R}^{n\times d}$ is the noise matrix, $U=\begin{bmatrix}\mathbb{u}_{1}&...&\mathbb{u}_{r}\end{bmatrix}\in\mathbb{R}^{
n\times r}$ has orthonormal column vectors being the latent $r$-dimensional factors or sample embeddings characterizing the underlying signal structure among the $n$ samples, and $Y=\begin{bmatrix}\sqrt{\lambda_{1}}\mathbf{y}_{1}&...&\sqrt{\lambda_{r}}
\mathbf{y}_{r}\end{bmatrix}\in\mathbb{R}^{d\times r}$ contains the feature loadings whose $(i,j)$ entry characterizing the weight of $j$th latent factor $\mathbb{u}_{j}$ in the $i$th feature. The above model essentially assumes that the data matrix $X$ contains a latent low-rank signal structure, which complies with many real applications and can be empirically verified by comparing the magnitude of the first few singular values with the other singular values. We assume the noise matrix $Z$ and the (rescaled) feature loading vectors $\{\mathbf{y}_{i}\}_{1\leq i\leq r}$ contain independent entries with zero mean and unit variance, but also remark that extensions to more general settings is possible (see discussions after Thm 1). Here, unlike classical theory of PCA where samples are assumed to be independent and features are correlated, we exchange the roles of the samples and features and model the underlying low-dimensional structure among samples by the latent factor $U$, or the low-rank correlation structure among the samples.

From the above model, the high-dimensional data matrix $X$ is a noisy realization of the low-dimensional latent signal matrix $U$ that encodes the true relationship among the samples. The true Euclidean angles between the noiseless samples $j$ and $k$ with respect to sample $i$ can be defined as
$\theta_{jk,i}=\arccos\big{(}\frac{(U_{j}-U_{i})\cdot(U_{k}-U_{i})}{\|U_{j}-U_{
i}\|\|U_{k}-U_{i}\|}\big{)}$,
where $U_{i}\in\mathbb{R}^{r}$ is the $i$th row of $U$, giving the true embedding of sample $i$. Our goal is to obtain reliable estimators of the latent Euclidean angles $\{\theta_{jk,i}\}$ based on the noisy data $X$. In our algorithm, we use the leading $r$ eigenvectors $\{\widehat{\mathbb{u}}_{i}\}_{1\leq i\leq r}$ of the Gram matrix $XX^{\top}$ (suppose the data is centered), and estimate $\theta_{jk,i}$ by
$\widehat{\theta}_{jk,i}=\arccos\big{(}\frac{(\widehat{U}_{j}-\widehat{U}_{i})
\cdot(\widehat{U}_{k}-\widehat{U}_{i})}{\|\widehat{U}_{j}-\widehat{U}_{i}\|\|
\widehat{U}_{k}-\widehat{U}_{i}\|}\big{)}$,
where $\widehat{U}_{i}\in\mathbb{R}^{r}$ is the $i$th row of $\widehat{U}=\begin{bmatrix}\widehat{\mathbb{u}}_{1}&...&\widehat{\mathbb{u}}_{
r}\end{bmatrix}$.
Our first result concerns the consistency of the latent angle estimation. For any pair $(i,j)$, we obtain the error bound for

$$ $|\widehat{U}_{i}^{\top}\widehat{U}_{j}-U_{i}^{\top}U_{j}|=|{\bf e}_{i}^{\top}( \widehat{U}\widehat{U}^{\top}-UU^{\top}){\bf e_{j}}|.$ (2) $$

The accuracy of estimating $U_{i}^{\top}U_{j}$ using $\widehat{U}_{i}^{\top}\widehat{U}_{j}$ is fundamental here since by

$$ $\frac{(U_{j}-U_{i})^{\top}(U_{k}-U_{i})}{\|U_{j}-U_{i}\|\|U_{k}-U_{i}\|}=\frac {U_{j}^{\top}U_{k}-U_{i}^{\top}U_{k}-U_{j}^{\top}U_{i}+\|U_{i}\|^{2}}{\sqrt{\| U_{j}\|^{2}+\|U_{i}\|^{2}-2U_{j}^{\top}U_{i}}\sqrt{\|U_{k}\|^{2}+\|U_{i}\|^{2} -2U_{k}^{\top}U_{i}}}$ (3) $$

the pairwise inner products $\{U_{a}^{\top}U_{b}:a,b\in\{i,j,k\}\}$ are the building blocks for the angle $\theta_{jk,i}$. In other words, the consistency of $\{U_{a}^{\top}U_{b}:a,b\in\{i,j,k\}\}$ implies the consistency of $\widehat{\theta}_{jk,i}$. Below we obtain the high-probability limit for the estimation error (2), which guarantees the estimation accuracy of $\widehat{\theta}_{jk,i}$ under sufficiently large signal-to-noise ratio.

To better present our results, we denote the aspect ratio $\phi=\frac{n}{d}$ and assume that $n^{1/C}\leq d\leq n^{C}$ for some constant $C>0$, characterizing the high dimensionality of the data. We define the rescaled Gram matrix $Q=\frac{1}{\sqrt{dn}}XX^{\top}$ and denote the population covariance $\Sigma=d^{-1}\mathbb{E}(XX^{\top})=I_{n}+UDU^{\top}=I_{n}+\phi^{1/2}\sum_{i=1}
^{r}\sigma_{i}\mathbb{u}_{i}\mathbb{u}_{i}^{\top}$, where $D=\text{diag}(\phi^{1/2}\sigma_{1},...,\phi^{1/2}\sigma_{r})$, and $\sigma_{1}\geq\sigma_{2}\geq...\geq\sigma_{r}>0$, so that $\{1+\phi^{1/2}\sigma_{i}\}_{1\leq i\leq r}$ are the leading $r$ eigenvalues of $\Sigma$.

###### Theorem 1 (Guarantee of spectral angle estimators) .

Suppose that $\sigma_{1}\geq...\geq\sigma_{r}\geq 1+c$ for some constant $c>0$.
Then for any $i,j\in\{1,2,...,n\}$, we have that, as $(n,d)\to\infty$,

$$ $|\widehat{U}_{i}^{\top}\widehat{U}_{j}-U_{i}^{\top}U_{j}|=\sum_{k=1}^{r}\frac{ (1+\phi^{1/2}\sigma_{k})u_{ki}u_{kj}}{\sigma_{k}(\sigma_{k}+\phi^{1/2})}+O_{P} (n^{-1/2+\epsilon}),$ (4) $$

for any small constant $\epsilon>0$, where we denote $\mathbb{u}_{k}=(u_{k1},...,u_{kn})$, for $1\leq k\leq r$.

From the above theorem as well as the relationship (3), we can see that the spectral angle estimator $\widehat{\theta}_{jk,i}$ can be arbitrarily close to $\theta_{jk,i}$ as the overall signal strength of the low-dimensional structure, characterized by the parameters $\{\sigma_{1},...,\sigma_{r}\}$, increases. Our analysis holds for general $\phi$, which may depend on $n$ and needs not to converge in $(0,\infty)$. In particular, our result implies the consistency of
$\widehat{\theta}_{jk,i}$ for any low-dimensional structures contained in $U$, that is, for any $\epsilon>0$, there exist sufficiently large $(\sigma_{1},...,\sigma_{r})$ such that

$$ $\lim_{n\to\infty}P(|\widehat{\theta}_{jk,i}-\theta_{jk.i}|>\epsilon)=0.$ (5) $$

We remark that the homoscedasticity assumption on the entries of the noise matrix $Z$ and the feature loading vectors $\{\mathbf{y}_{i}\}_{1\leq i\leq r}$ may be relaxed to more general settings, following the universality arguments in random matrix theory [].

Our next result concerns the non-negligible effect of high-dimensionality on the latent angle estimation. Here we assume $\phi$ remains bounded away from zero, that is, $\phi>c$ for some absolute constant $c>0$. We show that the angles between the original high-dimensional data points, that is,
$\bar{\theta}_{jk,i}:=\arccos\big{(}\frac{(X_{j}-X_{i})^{\top}(X_{k}-X_{i})}{\|
X_{j}-X_{i}\|\|X_{k}-X_{i}\|}\big{)}$,
can be substantially biased with respect to the latent angles.

###### Theorem 2 (Limitation of naive angle estimators) .

Under the assumption of Theorem 1, if we denote $X_{i}\in\mathbb{R}^{p}$ as the $i$th row of $X\in\mathbb{R}^{n\times d}$, it then holds that, for all $C>0$, there exist some $\Sigma$ with $\sigma_{1},...,\sigma_{r}>C$ such that, for any
distinct $i,j,k\in\{1,2,...,n\}$,
$\lim_{n\to\infty}P\left(\left|\bar{\theta}_{jk,i}-\theta_{jk,i}\right|\geq
\delta\right)=1,$
for some fixed constant $\delta>0$ that only depends on $\theta_{jk,i}$.

Comparing this with Eq. 5, we can see that for high-dimensional data, the naive angle estimators $\bar{\theta}_{jk,i}$ based on the original noisy high-dimensional observations can be substantially biased, regardless of signal strength. Theorems 1 and 2 together provide a theoretical justification and explain the
practical advantages of our spectral angle estimators for dealing with noisy high-dimensional data.

###### Theorem 1 (Guarantee of spectral angle estimators) .

###### Theorem 2 (Limitation of naive angle estimators) .

## 4 Experiments

To evaluate Mercat, we consider low-dimensional manifold datasets exemplifying unsolved issues in current LDEs, synthetic high-dimensional clustered data common in the literature, and real world applications including biology.
We compare against the state-of-the-art UMAP [], tSNE [], DensMAP [], and NCVis [].
To evaluate fairly, we consider a diverse set of metrics that measure how well properties of the high-dimensional data are preserved in the low-dimensional space.
We consider (i) distance preservation ($||.||$) measuring on reconstruction of all distances, (ii) preservation of angles ($\angle$) measuring how well the angles between data-points are preserved, (iii) neighborhood preservation ($\therefore$) measuring how well the closeness of local neighbors is preserved, and (iv) density preservation ($\odot$) measuring how well local sample density is preserved. All metrics are in the range $[-1,1]$, higher is better, and we provide all details and definitions in App. C.1.

For existing work, which focuses on reconstructing local distances correctly, the neighborhood parameter is crucial for embedding quality. We investigated the impact of this parameter on each method considering above metrics and found that there is a clear trade-off between local and global reconstruction, with increasing neighborhood respectively perplexity showing better reconstruction of long-range relationships but much worse reconstruction of local features, also evident in the embeddings. Most notably, there was no hyperparameter setting that was consistently better across metrics – different neighborhood sizes are optimal for different metrics.
We thus decided to stick to the recommended default if it was best for at least one metric, providing the analysis in App. C.4.

Figure: (a) Original data
Refer to caption: extracted/5666939/figs/smiley.png

Low-dimensional data
We first consider three datasets of 2 or 3 dimensions, as these can be directly visualized and hence compared to (see Supp. Sec. C.2 for details).
They exemplify difficult issues of current methodology for low-dimensional embeddings.
On simple data resembling a smiley face, a focus on reconstructing local distances results in the relative orientation of structures—here the eye, mouth and face outline—not being faithfully reconstructed and manifolds being distorted (see Fig. 2 top).
On a real 3D manifold representing the reconstruction of a mammoth [, ], we investigate how well complex manifolds are preserved, observing that current methods have issues deriving a meaningful embedding of the original data (see Fig. 2 mid row).
We observe the often discussed forming of ”arbitrary” clusters in UMAP, tSNE, and NCVis at varying degrees of intensity.
Studying different settings of the neighborhood parameter for existing work shows a trade-off between local and global feature reconstruction; with small neighborhood we see more clustering but also better capturing of local features, larger neighborhoods give better global reconstruction at the cost of fine-grained features (see App. C.4).
Mercat produces a faithful embedding of the mammoth that captures not only the main features but the pose of the animal.
For a simple circle in 2D, we observe that for such symmetric data current methods tend to break these symmetries deforming the circle and breaking it into clusters (see Fig. 2 bottom), which is consistent with the literature [].
We report all quantitative results in App. Tab. 3, observing that local neighborhoods are well preserved for existing methods, yet neither (global) distances nor angles are properly modeled. Both for a complex manifold as well as the highly symmetric circle data, the common objectives focusing on preserving local distances fail to yield faithful embeddings.
Mercat, on the other hand, yields embeddings that are as locally accurate as concurrent work, but outperforms them regarding distance and angle preservation (see Fig. 2, App. Tab. 3).

Cluster data
To evaluate on synthetic data that is standard in the literature, such as Gaussian mixtures, we consider five different datasets, varying number of clusters, distribution type, number of sample and density per cluster (see App. Sec. C).
We sample each dataset three times and report mean and standard deviation across different metrics in App. Tab. 3, giving visualizations for a fixed random seed in Supp. Fig. 6.
We observe that, consistent with the literature [], neither UMAP nor tSNE consistently outperform the other.
As expected, DensMAP, which explicitly optimizes for recovering local densities, outperforms all other methods on the investigated data in terms of density preservation.
Also, the general trends comparing between datasets are similar for all methods; the uniform data is more challenging than the simple Gaussian data (Unif5 vs Gauss5), and more clusters are harder to reconstruct (Gauss10 vs Gauss5).
On the challenging Gauss5-S and Gauss5-D data which have strongly varying densities between clusters, Mercat shows to be more robust than both tSNE and UMAP.
More importantly, across all experiments, we see that Mercat not only outperforms all competitors in terms of angle preservation—which it was optimized for—but also overall distance reconstruction and, perhaps surprisingly, preservation of density in most cases.
Overall, in all but one case our approach ranks first or second.

**Table 1: Real data results. We report angle preservation ($\angle$), distance preservation ($||.||$), neighborhood preservation ($\therefore$), and density preservation ($\odot$) between computed low-dimensional embeddings and original data on real data. All numbers are rounded to two decimal places, higher is better, and best method is in bold, second best is underlined.**
| Data | Metric | Mercat (ours) | UMAP | tSNE | NCVis | DensMAP |
| --- | --- | --- | --- | --- | --- | --- |
| Tab. Sap.<br>blood | $\angle$ | $\mathbf{.26}$ | $\underline{.17}$ | $.07$ | $.07$ | $.13$ |
| $||.||$ | $\mathbf{.25}$ | $\underline{.09}$ | $-.07$ | $-.01$ | $\underline{.09}$ |  |
| $\therefore$ | $\underline{.07}$ | $.02$ | $\mathbf{.29}$ | $.02$ | $.02$ |  |
| $\odot$ | $\mathbf{.22}$ | $.00$ | $.11$ | $-.16$ | $\underline{.21}$ |  |
| Murine<br>Pancreas<br>$n=50$ | $\angle$ | $\mathbf{.48}$ | $.33$ | $\underline{.46}$ | $.34$ | $.40$ |
| $||.||$ | $\mathbf{.61}$ | $.41$ | $.45$ | $.32$ | $\underline{.48}$ |  |
| $\therefore$ | $\underline{.10}$ | $.07$ | $\mathbf{.34}$ | $.06$ | $\underline{.10}$ |  |
| $\odot$ | $\mathbf{.56}$ | $-.22$ | $.01$ | $.14$ | $\underline{.27}$ |  |
| Hematop.<br>Paul et al.<br>$n=50$ | $\angle$ | $\mathbf{.86}$ | $.76$ | $\underline{.82}$ | $.34$ | $.77$ |
| $||.||$ | $\mathbf{.92}$ | $.75$ | $.81$ | $.44$ | $\underline{.82}$ |  |
| $\therefore$ | $\underline{.31}$ | $.28$ | $\mathbf{.35}$ | $.30$ | $.28$ |  |
| $\odot$ | $\mathbf{.66}$ | $.29$ | $.08$ | $.11$ | $\underline{.63}$ |  |
| MNIST<br>even<br>$n=50$ | $\angle$ | $\mathbf{.53}$ | $\underline{.35}$ | $.34$ | $.33$ | $.35$ |
| $||.||$ | $\mathbf{.61}$ | $.35$ | $\underline{.36}$ | $.34$ | $\underline{.36}$ |  |
| $\therefore$ | $.04$ | $\underline{.11}$ | $\mathbf{.20}$ | $\underline{.11}$ | $.10$ |  |
| $\odot$ | $.09$ | $-.06$ | $\underline{.14}$ | $.10$ | $\mathbf{.45}$ |  |
| Cell<br>Cycle<br>$n=50$ | $\angle$ | $\mathbf{.44}$ | $.20$ | $.24$ | $\underline{.28}$ | $.21$ |
| $||.||$ | $\mathbf{.51}$ | $.21$ | $.20$ | $\underline{.27}$ | $.21$ |  |
| $\therefore$ | $\underline{.18}$ | $.07$ | $\mathbf{.27}$ | $.07$ | $.11$ |  |
| $\odot$ | $.35$ | $\underline{.37}$ | $\mathbf{.43}$ | $.16$ | $.17$ |  |

Real world data
We evaluate the methods on three single-cell gene expression datasets of different origin resembling the most typical application of LDEs, in particular samples of human blood from the Tabula Sapiens project [], bone marrow in mice [], and from the Murine Pancreas []. We provide details on processing of the data in App. C.3. LDEs should capture the structure of blood cell differentiation.
We further consider the MNIST [], where we focus on even numbers, as state-of-the-art methods are presumably good at clustering and should hence be able to capture these well-separated classes better. Lastly, we consider a dataset of cells with estimated cell cycle stage [], an LDE can hence reflect the cyclic dependency of cell states.
We report results in Tab. 1 and provide all visualizations in App. Sec. C.7.

Consistent with our previous findings, we see that Mercat performs best in terms of angle and overall distance preservation. tSNE is best in reconstructing local neighborhoods, with Mercat usually taking second place.
As expected, DensMAP outperforms existing work in terms of reconstruction of (local) densities in most cases, with cell cycle data being an exception. Perhaps surprisingly, Mercat performs best in three out of five datasets regarding density reconstruction, despite not explicitly modeling this property.
In terms of quantitative results, Mercat seems to strike a balance of reconstructing both local as well as global structures also on real-world data.

Looking more closely into the Tabula Sapiens data (cf App. C.7), UMAP and DensMAP struggle with a proper fine-grained reflection of the data, as immune vs non-immune cells are dominating the overall structure and little structure is visible within immune cell clusters. tSNE learns several clusters, but dependencies between cell types are hard to make out. NCVis is able to find a more global structure, as well as differentiating locally between particular cell types, the visible global dependency looks, however, overly complex, much like the induced arbitrary bends on the Circle toy example (cf. Fig. 2d). Mercat learns a clearly visible and interpretable local and global structure reflecting relationships of different blood cell types, which together with the quantitative results indicate a more faithful reconstruction of the high-dimensional data.
On MNIST, we see the known exaggeration of clustering by existing methods, which gives a clearer separation of digits. Mercat shows a greater mixture of cluster boundaries. While this sacrifices a bit of local reconstruction, it seems to better represent global relationships (cf. Tab. 1). For this particular dataset, we observe a strong trade-off between local and global structure preservation.
Murine Pancreas as well as the human bone marrow data on a first glance look similar across methods, with all being able to distinguish cell types, encoding global dependencies that reflect hematopoiesis. Yet, tSNE and NCVis seem to have issues getting the long-range dependencies right, and all existing methods often show formations of seemingly arbitrary clusters.
On cell cycle data, only Mercat and DensMAP are able to capture the cyclic structure of the data, correctly embedding the dependencies between the different cell cycle stages. All other methods are not reflecting the cell stage transitions properly.

## 5 Discussion and Conclusion

In this work, we suggested a new paradigm for the computation of low-dimensional embeddings, arguing for a simpler approach compared to current methodologies. The central question we ask is whether reconstructing primarily local features, as common in state-of-the-art, is what we want, given that this approach profoundly constrains the quality of reconstruction of the global properties of the data.
Different from existing work, we cast the underlying optimization problem in terms of reconstructing angles between any set of three points correctly on a 2-dimensional sphere.
We suggested an efficient approach called Mercat that can easily learn LDEs by off-the-shelf gradient descent optimizers.
Further, we both empirically as well as theoretically motivate a sub-sampling approach and an initial denoising step, which improves the efficiency and robustness of the proposed algorithm for large and high-dimensional datasets. On synthetic, real-world, and easy-to-understand low-dimensional data, we show that our approach effectively recovers both local as well as global structures, outperforming existing methods despite, or maybe because of, its simplicity.
It thus supports the hypothesis that the trade-offs between local and global reconstruction are caused by algorithm choice rather than theoretical limitation.

While giving highly encouraging results, our work also leaves room for future improvements.
One direction of research could be further improvements of embedding quality; Mercat mostly outperforms existing work in terms of angle-, distance-, and neighborhood-preservation, yet is often seconded by DensMAP in terms of density preservation. While this may come by little surprise, as DensMAP explicitly optimizes for density preservation, it would still make for exciting future work to improve Mercat in that regard.
Also, algorithmic advances targeting the efficiency could be interesting; the current methodology of Mercat is applicable to arbitrary sized datasets as it linearly scales with the number of samples thanks to the subsampling procedure, but is not ideal due to a large constant factor. For close to online performance on very large datasets, similar to NCVis [] or FasttSNE [], additional work is required.
Lastly, we anticipate further theoretical insights, as the simple optimization loss lends itself for rigorous analysis.

## Appendix A Algorithm

### A.1 Computing geodesic angles with linear algebra

To efficiently compute geodesics on a sphere that is numerically stable and suitable for computation on graphics cards, we use the following trick.

Figure: Figure 3: Computing sphere angles with linear algebra. We visualize the idea of computing the angle $\alpha$ between two (geodesic) paths $\overline{AB},\overline{AC}$ on a sphere. The key insight is that the angle between the two geodesics is the same as the angle between the normals (visualized as arrows) of the two triangles $\Delta OAB$, $\Delta OAC$ in the ambient 3D space, with $O$ as center of the sphere.
Refer to caption: extracted/5666939/figs/GeodesicAngleTrick.png

To compute an angle $\angle BAC$ between the (geodesic) paths $\overline{AB}$, $\overline{AC}$ at point $A$, respecting the curvature of the sphere, we use the fact that the angle between these geodesics is the angle between the two planes $p_{OAB}$ and $p_{OAC}$ in the ambient 3D space, where $p_{ijk}$ is the plane that is spanned by the three points $i,j,k$, and $O$ is the center of the sphere, which we assume to be the origin of the space w.l.o.g..
Using the further insight that the angle between these two planes is the angle between their normal vectors, we can use the cross product to compute the two normal vectors, normalize the vectors to unit length and then compute the enclosed angle by using the definition of the scalar product

$$ $\displaystyle\angle BAC=\cos^{-1}\left(\frac{A\otimes B}{||A\otimes B||}\cdot \frac{A\otimes C}{||A\otimes C||}\right),$ $$

with $\otimes$ as the cross product.
We provide a visualization of this idea in Fig. 3.
In practice, as discussed in the main paper, we will drop the inverse cosine function in both high- and low-dimensional angle computations, which is a strictly monotone transformation.

## Appendix B Theory

Figure: Figure 4: Spectral analysis of angle space. For 500 samples randomly taken from human hematopoiesis data [] we show (a) the singular values of the matrix $\Theta_{i}$ of cosine-angles at sample $i$ (one line per sample) and (b) the distribution of effective rank of all $\Theta_{i}$ on this dataset. Angle matrices are of low (effective) rank, thus encourage subsampling of angles.
Refer to caption: extracted/5666939/figs/Rank_analysis.png

### B.1 Proof of Theorem 1

We define

$$ $f(\sigma_{i})=\frac{(1+\phi^{1/2}\sigma_{i})(1-\sigma_{i}^{-2})}{\phi^{1/2} \tau(\sigma_{i})}$ (6) $$

where $\tau(x)=\phi^{1/2}+\phi^{-1/2}+x+x^{-1}$.
By definition, we have

$$ $\widehat{U}_{i}^{\top}\widehat{U}_{j}={\bf e}^{\top}_{i}\widehat{U}\widehat{U} ^{\top}{\bf e}_{j},\qquad U_{i}^{\top}U_{j}={\bf e}^{\top}_{i}UU^{\top}{\bf e} _{j}.$ (7) $$

In the following lemma, proved in [, Section 5] and [, Section 5], concerns the limiting behavior of the bilinear form $\mathbb{w}_{1}^{\top}\widehat{U}\widehat{U}^{\top}\mathbb{w}_{2}$ for any unit vectors $\mathbb{w}_{1},\mathbb{w}_{2}\in\mathbb{R}^{n}$.

###### Lemma 1 .

Under the assumption of Theorem 1,
for any unit vectors $\mathbb{w}_{1},\mathbb{w}_{2}\in\mathbb{R}^{n}$, it holds that

$$ $\mathbb{w}_{1}^{\top}\widehat{U}\widehat{U}^{\top}\mathbb{w}_{2}=\sum_{k=1}^{r }f(\sigma_{k})\mathbb{w}_{1}^{\top}\mathbb{u}_{k}\mathbb{u}_{k}^{\top}\mathbb{ w}_{2}+O_{P}(n^{-1/2+\epsilon}),$ (8) $$

for any small constant $\epsilon>0$, where
$f(\sigma_{k})$ is defined in (6).

As a result, if we denote

$$ $u_{ki}={\bf e}_{i}^{\top}\mathbb{u}_{k},\qquad 1\leq i\leq n,\quad 1\leq k\leq r,$ $$

it then follows that

$$ $\displaystyle{\bf e}^{\top}_{i}\widehat{U}\widehat{U}^{\top}{\bf e}_{j}$ $\displaystyle=\sum_{k=1}^{r}u_{ki}u_{kj}f(\sigma_{k})+O_{P}(n^{-1/2+\epsilon}) ={\bf e}_{i}^{\top}U\Gamma U^{\top}{\bf e}_{j}+O_{P}(n^{-1/2+\epsilon})$ $$

where $\Gamma=\text{diag}\left(f(\sigma_{1}),...,f(\sigma_{r})\right)$.
As a result, it follows that

$$ $\displaystyle|{\bf e}^{\top}_{i}\widehat{U}\widehat{U}^{\top}{\bf e}_{j}-{\bf e }_{i}^{\top}UU^{\top}{\bf e}_{j}|$ $\displaystyle=|{\bf e}_{i}^{\top}U(\Gamma-I_{r})U^{\top}{\bf e}_{j}|+O_{P}(n^{ -1/2+\epsilon})$ $\displaystyle=\bigg{|}\sum_{k=1}^{r}u_{ki}u_{kj}\bigg{(}\frac{(1+\phi^{1/2} \sigma_{k})(1-\sigma_{k}^{-2})}{\phi^{1/2}\tau(\sigma_{k})}-1\bigg{)}\bigg{|}+ O_{P}(n^{-1/2+\epsilon})$ $\displaystyle=\sum_{k=1}^{r}\frac{(1+\phi^{1/2}\sigma_{k})u_{ki}u_{kj}}{\sigma _{k}(\sigma_{k}+\phi^{1/2})}+O_{P}(n^{-1/2+\epsilon})$ $$

This completes the proof.

###### Lemma 1 .

### B.2 Proof of Theorem 2

Note that $\Sigma=I+\phi\sum_{s=1}^{r}\sigma_{s}\mathbb{u}_{s}\mathbb{u}_{s}^{\top}$ implies

$$ $\Sigma_{ij}=\phi\sum_{s=1}^{r}\sigma_{s}u_{si}u_{sj}+\delta_{ij}.$ $$

Then we have

$$ $\displaystyle\frac{\Sigma_{jk}-\Sigma_{ik}-\Sigma_{ji}+\Sigma_{ii}}{\sqrt{ \Sigma_{jj}+\Sigma_{ii}-2\Sigma_{ji}}\sqrt{\Sigma_{kk}+\Sigma_{ii}-2\Sigma_{ki }}}$ $\displaystyle=\frac{\phi\sum_{s=1}^{r}\sigma_{s}(u_{sj}u_{sk}-u_{si}u_{sk}-u_{ si}u_{sj}+u_{si}^{2})+1}{\phi\sqrt{\sum_{s=1}^{r}\sigma_{s}(u_{i}-u_{j})^{2}} \sqrt{\sum_{s=1}^{r}\sigma_{s}(u_{i}-u_{k})^{2}}}$ $\displaystyle=\frac{(U_{j}-U_{i})^{\top}W(U_{k}-U_{i})+(\phi\sigma_{r})^{-1}}{ \sqrt{(U_{j}-U_{i})^{\top}W(U_{j}-U_{i})}\sqrt{(U_{k}-U_{i})^{\top}W(U_{k}-U_{ i})}},$ $$

where $W=\text{diag}(\sigma_{1}/\sigma_{r},\sigma_{2}/\sigma_{r}...,1)$. If we denote

$$ $\bm{\beta}=U_{j}-U_{i},\qquad\bm{\gamma}=U_{k}-U_{i},$ $$

and

$$ $\widetilde{\bm{\beta}}=W^{1/2}(U_{j}-U_{i}),\qquad\widetilde{\bm{\gamma}}=W^{1 /2}(U_{k}-U_{i}),$ $$

it follows that

$$ $\frac{\Sigma_{jk}-\Sigma_{ik}-\Sigma_{ji}+\Sigma_{ii}}{\sqrt{\Sigma_{jj}+ \Sigma_{ii}-2\Sigma_{ji}}\sqrt{\Sigma_{kk}+\Sigma_{ii}-2\Sigma_{ki}}}=\frac{ \widetilde{\bm{\beta}}^{\top}\widetilde{\bm{\gamma}}+(\sigma_{r}\phi)^{-1}}{\| \widetilde{\bm{\beta}}\|\|\widetilde{\bm{\gamma}}\|}$ (9) $$

On the other hand, we have

$$ $\displaystyle\cos\theta_{jk,i}$ $\displaystyle=\frac{(U_{j}-U_{i})^{\top}(U_{k}-U_{i})}{\|U_{j}-U_{i}\|\|U_{k}- U_{i}\|}$ $\displaystyle=\frac{\bm{\beta}^{\top}\bm{\gamma}}{\|\bm{\beta}\|\|\bm{\gamma} \|}.$ $$

Now if we denote $\theta=\angle(\bm{\beta},\bm{\gamma})$ and $\widetilde{\theta}=\angle(\widetilde{\bm{\beta}},\widetilde{\bm{\gamma}})$, it follows that

$$ $\displaystyle\bigg{|}\frac{\widetilde{\bm{\beta}}^{\top}\widetilde{\bm{\gamma} }+(\sigma_{r}\phi)^{-1}}{\|\widetilde{\bm{\beta}}\|\|\widetilde{\bm{\gamma}}\| }-\frac{\bm{\beta}^{\top}\bm{\gamma}}{\|\bm{\beta}\|\|\bm{\gamma}\|}\bigg{|}$ $\displaystyle\geq|\cos\widetilde{\theta}|-|\cos\theta|-\frac{1}{\phi\sigma_{r} \|\bm{\beta}\|\|\bm{\gamma}\|},$ $$

where in the last inequality we used $\|\bm{\beta}\|\leq\|\widetilde{\bm{\beta}}\|$ and $\|\bm{\gamma}\|\leq\|\widetilde{\bm{\gamma}}\|$.
To obtain the final result, we first show that, there exists some $W$ so that $\cos\widetilde{\theta}$ can be made arbitrarily close to $1$ or $-1$. Without loss of generality, we assume $\cos\theta>0$, and $\beta_{1}\gamma_{1}\neq 0$, where we use the notation $\bm{\beta}=(\beta_{1},...,\beta_{r})$ and $\bm{\gamma}=(\gamma_{1},...,\gamma_{r})$. Moreover, we denote $\alpha=\frac{1-\cos\theta}{\cos\theta}$ so that $1=(1+\alpha)\cos\theta$. Now if $\beta_{1}\gamma_{1}>0$, then we can always find $W$ so that $\sigma_{1}/\sigma_{r}$ is significantly larger than $\{\sigma_{2}/\sigma_{r},...,1\}$, and therefore either

$$ $\max\{\angle(\widetilde{\bm{\beta}},{\bf e}_{1}),\angle(\widetilde{\bm{\gamma} },{\bf e}_{1})\}<\frac{1}{2}\arccos\left(\bigg{(}1+\frac{\alpha}{2}\bigg{)} \cos\theta\right)$ $$

or

$$ $\max\{\angle(\widetilde{\bm{\beta}},-{\bf e}_{1}),\angle(\widetilde{\bm{\gamma }},-{\bf e}_{1})\}<\frac{1}{2}\arccos\left(\bigg{(}1+\frac{\alpha}{2}\bigg{)} \cos\theta\right)$ $$

holds. In either case, we have

$$ $\angle(\widetilde{\bm{\beta}},\widetilde{\bm{\gamma}})<\arccos\left(\bigg{(}1+ \frac{\alpha}{2}\bigg{)}\cos\theta\right)$ $$

so that

$$ $\cos\widetilde{\theta}>\bigg{(}1+\frac{\alpha}{2}\bigg{)}\cos\theta.$ $$

If instead $\beta_{1}\gamma_{1}<0$, then we can similarly choose $\sigma_{1}/\sigma_{r}$ sufficiently larger than $\{\sigma_{2}/\sigma_{r},...,1\}$ so that either

$$ $\max\{\angle(\widetilde{\bm{\beta}},{\bf e}_{1}),\angle(\widetilde{\bm{\gamma} },-{\bf e}_{1})\}<\frac{1}{2}\arccos\left(\bigg{(}1+\frac{\alpha}{2}\bigg{)} \cos\theta\right)$ $$

or

$$ $\max\{\angle(\widetilde{\bm{\beta}},-{\bf e}_{1}),\angle(\widetilde{\bm{\gamma }},{\bf e}_{1})\}<\frac{1}{2}\arccos\left(\bigg{(}1+\frac{\alpha}{2}\bigg{)} \cos\theta\right)$ $$

holds. In either case, we have

$$ $\angle(\widetilde{\bm{\beta}},\widetilde{\bm{\gamma}})>\pi-\arccos\left(\bigg{ (}1+\frac{\alpha}{2}\bigg{)}\cos\theta\right),$ $$

so that

$$ $\cos\widetilde{\theta}<-\bigg{(}1+\frac{\alpha}{2}\bigg{)}\cos\theta.$ $$

As a result, we have

$$ $\bigg{|}\frac{\widetilde{\bm{\beta}}^{\top}\widetilde{\bm{\gamma}}+(\sigma_{r} \phi)^{-1}}{\|\widetilde{\bm{\beta}}\|\|\widetilde{\bm{\gamma}}\|}-\frac{\bm{ \beta}^{\top}\bm{\gamma}}{\|\bm{\beta}\|\|\bm{\gamma}\|}\bigg{|}\geq\frac{ \alpha}{2}\cos\theta-\frac{1}{\phi\sigma_{r}\|\bm{\beta}\|\|\bm{\gamma}\|}.$ (10) $$

Finally, with $W$ and $U$ fixed and $\phi$ bounded away from zero, we can always choose sufficiently large $\sigma_{r}>0$ such that

$$ $\frac{1}{\phi\sigma_{r}\|\bm{\beta}\|\|\bm{\gamma}\|}<\frac{\alpha}{4}\cos\theta.$ $$

Combining the above results, we have

$$ $\bigg{|}\frac{\Sigma_{jk}-\Sigma_{ik}-\Sigma_{ji}+\Sigma_{ii}}{\sqrt{\Sigma_{ jj}+\Sigma_{ii}-2\Sigma_{ji}}\sqrt{\Sigma_{kk}+\Sigma_{ii}-2\Sigma_{ki}}}-\cos \theta_{jk,i}\bigg{|}\geq\frac{\alpha}{4}\cos\theta,$ (11) $$

or

$$ $\bigg{|}\arccos\left(\frac{\Sigma_{jk}-\Sigma_{ik}-\Sigma_{ji}+\Sigma_{ii}}{ \sqrt{\Sigma_{jj}+\Sigma_{ii}-2\Sigma_{ji}}\sqrt{\Sigma_{kk}+\Sigma_{ii}-2 \Sigma_{ki}}}\right)-\theta\bigg{|}\geq\delta,$ (12) $$

for some constant $\delta>0$ only depending on $\theta$.
Finally, it suffices to note that, by the law of large numbers and the continuous mapping theorem, we have

$$ $\arccos\bigg{(}\frac{(X_{j}-X_{i})\cdot(X_{k}-X_{i})}{\|X_{j}-X_{i}\|\|X_{k}-X _{i}\|}\bigg{)}\to_{P}\arccos\left(\frac{\Sigma_{jk}-\Sigma_{ik}-\Sigma_{ji}+ \Sigma_{ii}}{\sqrt{\Sigma_{jj}+\Sigma_{ii}-2\Sigma_{ji}}\sqrt{\Sigma_{kk}+ \Sigma_{ii}-2\Sigma_{ki}}}\right).$ (13) $$

This along with (12) completes the proof of the theorem.

### B.3 Efficacy of subsampling

Here, we provide some theoretical insights that partially explains the efficacy of our subsampling procedure. Recall that at each optimization iteration, for each data point $i$, instead of using of all the entries in the angle matrix $M_{i}=(\angle jik)_{1\leq j,k\leq n}$, we only take a random subset of the entries. Our hope is that such a random subset contains sufficient information about the whole matrix. This is in the same spirit as the matrix completion problem where the goal is to recover the missing matrix entries from a small number of randomly observed entries [, , , ]. From theory of matrix completion, a critical condition enabling precise local-to-global reconstruction is known as the incoherence condition, which essentially requires that the matrix is approximately low-rank and its leading singular vectors are relatively “spread out,” effectively avoiding any outliers in the data matrix.
In our case, the spiked population model automatically implies the approximate low-rankness of the cosine-angle matrix $\widehat{\Theta}_{i}=(\widehat{\theta}_{jk,i})_{1\leq j\neq k\leq n}$, which follows from (3) and that

$$ $\Theta_{i}\equiv(\theta_{jk,i})_{1\leq j,k\leq n}=\bigg{(}\frac{(U_{j}-U_{i})^ {\top}(U_{k}-U_{i})}{\|U_{j}-U_{i}\|\|U_{k}-U_{i}\|}\bigg{)}=D^{-1/2}VV^{\top} D^{-1/2},$ (14) $$

where

$$ $V=\begin{bmatrix}(U_{1}-U_{i})^{\top}\\ (U_{2}-U_{i})^{\top}\\ ...\\ (U_{n}-U_{i})^{\top}\end{bmatrix}\in\mathbb{R}^{n\times r},\quad D=\text{diag} (\|U_{1}-U_{i}\|^{2},...,\|U_{n}-U_{i}\|^{2}),$ (15) $$

showing that $\Theta_{i}$ has rank at most $r$. If we denote $W\in\mathbb{R}^{n\times r}$ as the matrix of singular vectors of $\Theta_{i}$, the incoherence condition amounts to saying that

$$ $\bigg{\|}WW^{\top}-\frac{r}{n}I_{n}\bigg{\|}_{\max}\leq\mu\frac{\sqrt{r}}{n}$ (16) $$

for some small constant $\mu>0$, where $\|(a_{ij})\|_{\max}=\max_{i,j}|a_{ij}|$. In particular, the incoherence condition (16) is likely satisfied if the low-dimensional signal structure with respect to the $i$th data point, encoded by $\{U_{j}-U_{i}\}_{1\leq j\leq n}$, has certain smoothness property and does not contain outliers deviating significantly from the bulk, which is the case for many applications. For example, in typical biological applications an outlier removal is part of the preprocessing pipeline [, ].

## Appendix C Experiments

### C.1 Computation of evaluation metrics

In the following, we provide an overview on how the evaluation metrics are defined.

- •
distance preservation ($||.||$) measured as Spearman Rank correlation coefficient between high- and low-dimensional distances, capturing how well overall structure is preserved. Distances for Mercat embeddings are computed from geodesics on the sphere.
- •
neighborhood preservation ($\therefore$) as measured by the mean jaccard index of the $k$-nearest neighbors (here, $k=50$) in high- and low-dimensional space across all points, $1/n\sum_{i}\frac{|knn(X,i)\cap knn(Y,i)|}{|knn(X,i)\cup knn(Y,i)|}$, where $knn(X,i)$ gives the indices of the $k$ nearest neighbors in $X$, capturing how accurate local structures are embedded. Before neighborhood computation, we denoise using ScreeNOT [].
- •
density preservation ($\odot$), which reflects how well differences in densities are captured in the embedding, a recent point of interest in the literature [, ]. We measure this by comparing the number of points that fall in spheres of constant radius around each point. More concretely, we compute the average distance of the $25th$-nearest neighbor in high- and low-dimensional space, $\bar{k}_{high}$ and $\bar{k}_{low}$, and for each sample $i$ compute the local density as number of points that fall into a sphere centered at $i$ of radius $\bar{k}_{high}$ resp. $\bar{k}_{low}$. The Pearson correlation coefficient between the obtained sphere densities gives our final metric.
- •
preservation of angles ($\angle$) between any three points measured as the Pearson Correlation coefficient between angles in high- and low-dimensional space, which captures how well global relationships, such as orientation of clusters are preserved. For practical purposes, as this computation is cubic in the number of points, we again sample for each point $i$ $64$ other points at random and compute the angle at $i$ and all combination of other points.

### C.2 Reproducibility – Generation of data

#### Smiley

To obtain the Smiley dataset, we sample $n=3000$ points as follows.
A quarter of these points are used for the eyes, where we first draw a radius for each point as $e^{\prime}_{r}\sim U(0,1)$ and further transform this radius to get $e_{r}=.1\sqrt{e^{\prime}_{r}}$. We additionally draw an angle $e_{\theta}\sim U(0,2\pi)$.
The actual points are then assigned to the 2D coordinates $x=e_{r}\sin(e_{\theta}),~{}y=e_{r}*\cos(e_{\theta})$. Half of these samples are then offset by $(.25,.25)$, the other half by $(-.25,.25)$, resulting in the final coordinates of the eyes.
For the face outline we dedicate half of the overall points, first sampling a radius $f^{\prime}_{r}\sim U(.9^{2},1)$, which is transformed to get $f_{r}=\sqrt{f^{\prime}_{r}}$.
We further draw an angle $f_{\theta}\sim U(0,2\pi)$ and compute the final coordinates as
$x=f_{r}\sin(f_{\theta}),~{}y=f_{r}*\cos(f_{\theta})$.
Lastly, we dedicate the remaining (quarter of) points to the mouth, sampling $m^{\prime}_{r}\sim U(.45^{2},.55^{2})$, which is transformed to get $m_{r}=\sqrt{m^{\prime}_{r}}$.
We further draw an angle $m_{\theta}\sim U(0,\pi)$ and compute the 2D coordinates as $x=m_{r}\sin(m_{\theta}),~{}y=-m_{r}*\cos(m_{\theta})$.
Lastly, we scale the whole data by 2, concluding the data generation process

#### Circle

For the Circle data, we sample $n=900$ angle $c_{\theta}\sim U(0,2\pi)$ and compute the original circle as
$x=3\cos(c_{\theta}),~{}y=3*\sin(c_{\theta})$. We then add iid noise sampled from $N(0,.01)$ to both dimensions.

#### Generation of synthetic data

We generate (i) Unif5, a dataset in 50 dimensions of 5 uniform clusters with $100$ samples each, with each dimension iid from $U(0,1)$ and different centers sampled from $U(-10,10)$, (ii) Gauss5, a dataset in 50 dimensions from 5 Gaussians with mean $\mu$ sampled from $U(-10,10)$ (iid for each dimension) and standard deviation $\sigma$ sampled from $U(.5,2)$ (iid for each dimension, all dimensions have covariance of 0), each cluster having 100 samples each, (iii) Gauss10, a dataset in 50 dimensions from 10 Gaussians with mean $\mu$ sampled from $U(-10,10)$ (iid for each dimension) and standard deviation $\sigma$ sampled from $U(.5,2)$ (iid for each dimension, all dimensions have covariance of 0), each cluster having 100 samples each, (iv) Gauss5-S, which is generated similar as Gauss5, but with different number of samples per cluster, namely 50,100,150,200, and 250 samples, and (v) Gauss5-D, which is generated similar as Gauss5, but with different densities per cluster using a covariance matrix as a diagonal matrix where entries are set 1, 2, 3, 4, and 5 for each cluster respectively.

### C.3 Reproducibility – Preprocessing of real data

##### Tabula Sapiens human blood

We obtained the human blood samples from the Tabula Sapiens project through the CZ CELLxGENE portal, preprocessed as Seurat object. We proceeded by filtering for data from the 10x 3’ v3 assay to avoid strong batch effects due to different sequencing platforms.
To filter for protein-coding genes – excluding genes encoded in the mitochondrium – we used the Gencode v38 genome annotation.
We further filtered for genes that were expressed in at least one sample (i.e., sum of gene expression across samples was greater than zero).
The annotated cell type in the data object was used for labeling.

##### Murine pancreas

We obtained pre-processed single-cell gene expression data through the Gene Expression Omnibus (accession id GSE132188).
To filter for protein-coding genes, we used the genome annotation GRCm39.110.
As before, we further filtered for genes that were expressed in at least one sample.
For cell annotation, we use the provided clusters used in Figure 3 of the original publication [].

##### Mouse bone marrow

We obtain the pre-processed single-cell data of Paul et al. [] from the PAGA repository(^2^22https://github.com/theislab/paga) [].

##### Cell cycle data

The HeLa cell cycle annotated data was obtained following the github repository(^3^33https://github.com/danielschw188/Revelio) of the original authors [], using the estimated phase as labels.

### C.4 Hyperparameter choices

We checked different hyper-parameter settings for existing work, focusing on varying the neighborhood respectively perplexity scores for UMAP, tSNE, NCVis, and DensMAP, as this is known to be one of the most deciding factors of embedding quality [].
As datasets, we consider a representative subset using Unif5 from the cluster datasets, Mammoth from the low-dimensional manifold datasets, for both of which we vary the parameter $\theta\in{15,30,50,100,200}$, and hematopoiesis data of Paul et al. from the real world datasets, for which we consider $\theta\in{15,30,100,200,500}$, as it is considerably larger.
We give the quantitative results in Tab. 2 and provide visualization of the mammoth reconstructions in Fig. 5, as we can compare these with the visualization of the original data (cf Fig. 2).

Across data, we see that quantitatively there is no single best parameter $\theta$, not across datasets, but more importantly, not within a dataset: varying the locality parameter $\theta$ (neighborhood or perplexity) means trading off local reconstruction performance against global reconstruction performance.
This also becomes evident in the visualizations for mammoth (Fig. 5), where for UMAP and DensMAP, which arguably give better reconstructions than competing methods, at smaller neighborhood size parameters the shape of the hip or leg bones as well as ribcage are still visible, at higher resolution the overall global structure looks like a more natural animal pose (albeit still wrong).
We, hence, decided to use the recommended default neighborhood parameter if at least one metric was ”optimal” during our evaluation. All other parameters were kept at their default value, noting that training converged in all but one case.
This particular case was UMAP on the Tabula Sapiens blood data, where training with the default parameter yielded a particularly bad, artifacted visualization (albeit decent performance on local reconstruction). We then decided to set the neighborhood parameter to $50$ to arrive at a meaningful embedding.
For all remaining experiments we use the following setting:

- UMAP
$\textbf{n_neighbor=15 (recommended default)};$ use spectral initialization; $min_dist=0.1;$
- tSNE
$\textbf{perplexity = 30 (recommended default)};$ $initial_dims=50;$ $theta=0.5;$ use PCA initialization;$max_iter=1000;$ normalize data; $momentum=0.5;$ $final_momentum=0.8;$ $eta=200;$ $exaggeration_factor=12$
- NCVis
n_neighbors=15 (recommended default); $n_epochs=50$; $n_init_epochs=20$; $min_dist=0.4$
- DMAP
n_neighbors = 30 (recommended default); spectral initialization, $dens_frac=0.3$; $dens_{l}ambda=0.1$; $dens_var_shift=0.1$; $n_epochs=750$; $learning_rate=1$; $min_dist=0.1$

For Mercat, we use the standard parameters for the Adam optimizer as recommended in the original paper []. Throughout all experiments we set the initial learning rate to $0.01$, and have a multiplicative learning rate schedule $\gamma$, multiplying by $0.1$ at predefined iterations (i.e., reducing the learning rate by an order of magnitude). As discussed in the main paper, we use an angle subsampling of $64$, and a batch size of $64$. For all synthetic and toy experiments, we run for $t=1000$ iterations, with a learning rate change at $\gamma=[350]$.

For real world data we reduce the number of iterations, as we do a batched learning approach and hence need much fewer iterations to see the same number of samples (and hence angles) as in the synthetic case studies. In particular, for MNIST we use $t=250,\gamma=[100]$, for Tabula Sapiens and Murine Pancreas we use $t=50,\gamma=[10,30]$, for human bone marrow and cell cycle data we use $t=200,\gamma=[50,150]$.

Note that in principle it is possible to optimize these hyperparameters (learning rate, batch size, subsampling, etc) to further improve Mercat embeddings by calibrating based on angle reconstruction. We instead wanted to keep parameters constant across experiments to show Mercat’s wide applicability with a standard set of parameters and only vary the number of iterations and learning rate schedule linked to these iterations.

Figure: Table 2: Results on different neighborhood parametrization. We report angle preservation ($\angle$), distance preservation ($||.||$), neighborhood preservation ($\therefore$), and density preservation ($\odot$) between computed low-dimensional embeddings and original data. All numbers are rounded to two decimal places, higher is better, and best result is in bold.

Figure: (a) UMAP 15
Refer to caption: extracted/5666939/figs/umap_mam_hyper15.png

### C.5 Visualization-optimal rotations for 2D conformal maps

For a Mercat embedding, for any rotation or translation on the sphere, the embeddings obviously are equal, both in terms of loss and any distance or angle-based metrics on the sphere.
However, for visualization on a 2D map, such as a Mercator projection, which is a conformal map of the sphere, points close to the equator of this map show much less distortion in terms of distances compared to points close to the pole.
This can be seen in for example maps of planet earth commonly used in an atlas or most other print media, where the arctis or antarctis appear extremely stretched—or overly large—compared to their actual size relative to e.g. Europe.
For 2D visualizations of any Mercat embedding $Y$, we hence use a rotation that puts as many points as possible close to the equator, thus avoiding as much ”stretching” as possible.
To this end, we compute a simple grid of rotation angles $\alpha\in[-\pi/2,\pi/2],\beta\in[0,\pi]$ with a granularity of 40 (i.e., grid values in steps of $\pi/40$) for rotation matrix $R_{\alpha,\beta}=R_{\alpha}R_{\beta}$, with

$$ $\displaystyle R_{\alpha}=\begin{pmatrix}\cos(\alpha)&0&-\sin(\alpha)\\ 0&1&0\\ \sin(\alpha)&0&\cos(\alpha)\\ \end{pmatrix},$ $\displaystyle R_{\beta}=\begin{pmatrix}\cos(\beta)&-\sin(\beta)&0\\ \sin(\beta)&\cos(\beta)&0\\ 0&0&1\end{pmatrix}.$ $$

By evaluating a simple penalty based on the sum of squared latitudes across all points in the rotated embedding $Y^{r}=YR_{\alpha,\beta}$,
defined as $\sum_{i}(|\cos^{-1}(Y^{r}_{i})-\pi/2|)^{2}$,
we can optimize for a equator-favoring rotation for visualization purposes.
We use this approach to generate any 2D maps of Mercat embeddings.

### C.6 Synthetic data results

We give visualizations of the generated embeddings for synthetic data in Fig. 6 and quantitative evaluation in Tab. 3. All visualizations are for seed $1$ of the repeated experiments, results are visually very similar across seeds, as also evident from the performance metrics.

**Table 3: Synthetic and toy data results. We report angle preservation ($\angle$), distance preservation ($||.||$), neighborhood preservation ($\therefore$), and density preservation ($\odot$) between computed low-dimensional embeddings and original data on synthetic benchmarks. We report mean and standard deviation across $3$ repetitions of data generation, except mammoth and smiley which are static data sets. All numbers are rounded to two decimal places, higher is better, and best method is in bold, second best is underlined.**
| Data | Metric | Mercat (ours) | UMAP | tSNE | NCVis | DensMAP |
| --- | --- | --- | --- | --- | --- | --- |
| Smiley | $\angle$ | $\mathbf{1.0}$ | $.10$ | $.25$ | $.14$ | $\underline{.26}$ |
| $||.||$ | $\mathbf{1.0}$ | $.11$ | $\underline{.39}$ | $.22$ | $.37$ |  |
| $\therefore$ | $\mathbf{.85}$ | $.79$ | $\underline{.84}$ | $.75$ | $.83$ |  |
| $\odot$ | $\mathbf{.98}$ | $-.16$ | $-.32$ | $.27$ | $\underline{.88}$ |  |
| Mammoth | $\angle$ | $\mathbf{.95}$ | $.56$ | $.50$ | $.16$ | $.68$ |
| $||.||$ | $\mathbf{.99}$ | $.75$ | $.61$ | $.21$ | $.88$ |  |
| $\therefore$ | $.31$ | $.57$ | $\mathbf{.65}$ | $.54$ | $\underline{.60}$ |  |
| $\odot$ | $\underline{.59}$ | $.01$ | $.10$ | $.31$ | $\mathbf{.73}$ |  |
| Circle | $\angle$ | $\mathbf{.99}$$\mathbf{\pm.00}$ | $.73$$\pm.00$ | $.64$$\pm.00$ | $.28$$\pm.03$ | $\underline{.95}$$\underline{\pm.01}$ |
| $||.||$ | $\mathbf{.99}$$\mathbf{\pm.00}$ | $.85$$\pm.00$ | $.72$$\pm.00$ | $.44$$\pm.06$ | $\underline{.96}$$\underline{\pm.01}$ |  |
| $\therefore$ | $\mathbf{.90}$$\mathbf{\pm.00}$ | $.83$$\pm.00$ | $\mathbf{.90}$$\mathbf{\pm.00}$ | $.77$$\pm.03$ | $\mathbf{.90}$$\mathbf{\pm.01}$ |  |
| $\odot$ | $\underline{.77}$$\underline{\pm.00}$ | $.10$$\pm.0$ | $.47$$\pm.00$ | $.20$$\pm.13$ | $\mathbf{.89}$$\mathbf{\pm.01}$ |  |
| Unif5 | $\angle$ | $\mathbf{.67}$$\mathbf{\pm.01}$ | $.49$$\pm.02$ | $.50$$\pm.04$ | $\underline{.51}$$\underline{\pm.02}$ | $\underline{.51}$$\underline{\pm.02}$ |
| $||.||$ | $\mathbf{.90}$$\mathbf{\pm.05}$ | $.41$$\pm.07$ | $.53$$\pm.20$ | $.44$$\pm.13$ | $\underline{.58}$$\underline{\pm.12}$ |  |
| $\therefore$ | $\mathbf{.49}$$\mathbf{\pm.02}$ | $.36$$\pm.01$ | $\underline{.37}$$\underline{\pm.01}$ | $\underline{.37}$$\underline{\pm.01}$ | $\underline{.37}$$\underline{\pm.00}$ |  |
| $\odot$ | $.22$$\pm.02$ | $.18$$\pm.03$ | $\underline{.59}$$\underline{\pm.03}$ | $.45$$\pm.03$ | $\mathbf{.61}$$\mathbf{\pm.04}$ |  |
| Gauss5 | $\angle$ | $\mathbf{.72}$$\mathbf{\pm.00}$ | $.53$$\pm.02$ | $.50$$\pm.00$ | $\underline{.57}$$\underline{\pm.01}$ | $.49$$\pm.04$ |
| $||.||$ | $\mathbf{.93}$$\mathbf{\pm.00}$ | $.66$$\pm.00$ | $.52$$\pm.00$ | $\underline{.69}$$\underline{\pm.01}$ | $.45$$\pm.12$ |  |
| $\therefore$ | $\mathbf{.49}$$\mathbf{\pm.00}$ | $.37$$\pm.00$ | $\underline{.38}$$\underline{\pm.00}$ | $\underline{.38}$$\underline{\pm.00}$ | $.37$$\pm.00$ |  |
| $\odot$ | $.33$$\pm.00$ | $.15$$\pm.00$ | $\underline{.59}$$\underline{\pm.00}$ | $.46$$\pm.04$ | $\mathbf{.65}$$\mathbf{\pm.00}$ |  |
| Gauss10 | $\angle$ | $\mathbf{.61}$$\mathbf{\pm.00}$ | $.33$$\pm.00$ | $\underline{.35}$$\underline{\pm.00}$ | $\underline{.35}$$\underline{\pm.00}$ | $\underline{.35}$$\underline{\pm.00}$ |
| $||.||$ | $\mathbf{.82}$$\mathbf{\pm.00}$ | $\underline{.26}$$\underline{\pm.00}$ | $.17$$\pm.00$ | $.21$$\pm.02$ | $.21$$\pm.08$ |  |
| $\therefore$ | $\mathbf{.44}$$\mathbf{\pm.00}$ | $.37$$\pm.00$ | $\underline{.40}$$\underline{\pm.00}$ | $\underline{.40}$$\underline{\pm.00}$ | $.38$$\pm.00$ |  |
| $\odot$ | $.22$$\pm.00$ | $.09$$\pm.00$ | $\underline{.62}$$\underline{\pm.00}$ | $.56$$\pm.00$ | $\mathbf{.70}$$\mathbf{\pm.01}$ |  |
| Gauss5-S | $\angle$ | $\mathbf{.70}$$\mathbf{\pm.00}$ | $\underline{.53}$$\underline{\pm.00}$ | $.46$$\pm.00$ | $.50$$\pm.00$ | $.51$$\pm.01$ |
| $||.||$ | $\mathbf{.90}$$\mathbf{\pm.00}$ | $\underline{.61}$$\underline{\pm.00}$ | $.38$$\pm.00$ | $.44$$\pm.00$ | $.52$$\pm.27$ |  |
| $\therefore$ | $\mathbf{.36}$$\mathbf{\pm.00}$ | $.24$$\pm.00$ | $\underline{.26}$$\underline{\pm.00}$ | $.25$$\pm.00$ | $.25$$\pm.00$ |  |
| $\odot$ | $\underline{.38}$$\underline{\pm.02}$ | $-.06$$\pm.0$ | $.23$$\pm.00$ | $.33$$\pm.05$ | $\mathbf{.57}$$\mathbf{\pm.01}$ |  |
| Gauss5-D | $\angle$ | $\mathbf{.69}$$\mathbf{\pm.00}$ | $.49$$\pm.00$ | $.49$$\pm.00$ | $\underline{.57}$$\underline{\pm.01}$ | $.50$$\pm.01$ |
| $||.||$ | $\mathbf{.88}$$\mathbf{\pm.00}$ | $.49$$\pm.00$ | $.56$$\pm.00$ | $\underline{.76}$$\underline{\pm.01}$ | $.59$$\pm.02$ |  |
| $\therefore$ | $\mathbf{.51}$$\mathbf{\pm.00}$ | $.36$$\pm.00$ | $\underline{.38}$$\underline{\pm.00}$ | $.37$$\pm.00$ | $.36$$\pm.01$ |  |
| $\odot$ | $\underline{.60}$$\underline{\pm.00}$ | $-.15$$\pm.0$ | $.06$$\pm.00$ | $.05$$\pm.03$ | $\mathbf{.74}$$\mathbf{\pm.01}$ |  |

Figure: (a) UMAP
Refer to caption: extracted/5666939/figs/XU_umap.png

### C.7 Visualizations for real world data

We provide visualizations of the embeddings generated by all methods on real data in Fig.  7, 8, 9, 10, 11 and give runtime estimates in Tab. 4, all methods being run on the same commodity hardware (CPU: 13th Gen. Intel Core i5-1350P, RAM: 32GB DDR5 5600MHz, OS: Debian 12).

Figure: (a) $\textsc{UMAP}^{*}$
Refer to caption: extracted/5666939/figs/blood_umap.png

Figure: (a) UMAP
Refer to caption: extracted/5666939/figs/murinepanc_umap.png

Figure: (a) UMAP
Refer to caption: extracted/5666939/figs/paul_umap.png

Figure: (a) UMAP
Refer to caption: extracted/5666939/figs/cc_umap.png

Figure: (a) UMAP
Refer to caption: extracted/5666939/figs/mnist_umap.png

**Table 4: Runtime on real data. We report wall clock running time for all methods. *We use a standard implementation of tSNE available in R and note that [] proposed a faster version of tSNE. A comparison of all runtime-improvements for standard LDE approaches is out of scope of this paper, as runtime efficiency is not the prime interest.**
| Data | Mercat (ours) | Mercat GPU | UMAP | tSNE * | NCVis | DensMAP |
| --- | --- | --- | --- | --- | --- | --- |
| Tab. Sap.<br>blood | $58.8m$ | $24.7m$ | $58.0m$ | $10.2m$ | $17s$ | $60s$ |
| Murine<br>Pancreas | $2.2h$ | $20.6m$ | $12s$ | $6s$ | $11s$ | $47s$ |
| Hematop.<br>Paul et al. | $17.7m$ | $18.8m$ | $18s$ | $12s$ | $1s$ | $28s$ |
| Cell Cycle | $9.1m$ | $6.5m$ | $6s$ | $3s$ | $1s$ | $16s$ |
| MNIST even | $4.2h$ | $3.5h$ | $1.8m$ | $1.7m$ | $2s$ | $46s$ |

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: The idea behind the paper is interesting, and the authors even try to convey theoretically how to extract relevant information from the given dataset (angles between triplets of points). However, from the main text it is hard to find clear results why angle preserving embedding is preferable over the current embedding techniques. I did not see a clear theorem or clear empirical results on real datasets that show the preferability of such embedding. In the next lines I will go over some of the problems I had with the main text-

1. General comment- From my understanding of the paper you are working on visualization problems which is more specific than the low dimensional embedding algorithms. 

2. The Related work is one big paragraph of 32 lines. It would be preferable if this would have been broken into multiple paragraphs.

3. Initialization- The proposed method is initialized using the projection of the first two PC components onto half a sphere far away from the poles. It would be beneficial if the authors can discuss why is PCA a good initialization with respect to angles. For example, if the data is a 2-dimensional swiss roll that is embedded in 3-dimensional space I am not sure that PCA is preferable.

4. Theoretical - It will be beneficial if the authors will add a theorem that describes which data structures can be embedded based on angles.  For example, in manifold learning the intrinsic dimension plays a main role when considering to embed data.

5. Theoretical- It would be beneficial to a proof for for the conclusion that appears in Equation 4 (derived from equation 2 and theorem 2).

6. Experiments- Each evaluation metric is described in a few words and only in the Appendix they describe it technically. It will be beneficial to discuss them a bit.

7. Experiments- The proposed algorithm embeds the data in 3-dimensional space with some restrictions, while it compares to visualization techniques that use 2-dimensional space. Some of the compared algorithm can be used to generate a 3-dimensional embedding, it will be beneficial to compare to their 3-dimensional embedding version.

8. Experiments- It will be beneficial to add the PCA embedding in the different experiments, as the proposed method uses it as initialization. This will allow the reader to understand the adjustments made by your method to it. For example, in Fig 2. /"Low Dimensional data" experiment - the embedding of the proposed method seems to be close to the PCA embedding.

9. Experiments- I think that it will be beneficial to show the embedding of spectral algorithms such as Laplacian Eigenmaps and Diffusion maps.

10. Experiments-  Cluster data- The authors considered experiments on Cluster Data within the main text. However, the visualization and the tables with the performance measure were left out in the Appendix.  It will be beneficial if the authors add a figure/table about it in the main text.

11. Experiments- Real data- Some of these dataset have labels. It will be beneficial to see how well the different embeddings techniques embeds samples with similar labels together. I think that such visualization/ performance measures will be more beneficial than current performance measures.
soundness: 2
strengths: See questions.
summary: The authors propose a new technique that embeds a given dataset on  a 2-dimensional sphere within 3-dimensional space, while approximately preserving the angles between any triplets of points.
weaknesses: See questions.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: * PCA is used in Algorithm 1 to get \hat{X}. Much high-frequency information is lost in this step.
* Subsampling can lose much information about the manifold.
* Why do you multiply with 0.7*pi in the algorithm?
* Where are the definitions of performance metrics (angle preservation, distance preservation, neighborhood preservation, and density preservation)?
* A reference and discussion is missing. Faithfulness has already been discussed in the following paper.

J. Venna et al. Information Retrieval Perspective to Nonlinear Dimensionality Reduction for Data Visualization. In JMLR 2010.
soundness: 1
strengths: * An alternative way for nonlinear dimensionality reduction
* Some theoretical guarantees are provided
summary: The paper presents a low-dimensional embedding method named MERCAT that preserves angular relationships to maintain both local and global data structures. Unlike traditional methods that prioritize local distances and often distort global structures, MERCAT maps high-dimensional data onto a 2D unit sphere by optimizing angles between points. This approach is claimed to achieve consistent structure preservation across scales and advantages over methods like UMAP and tSNE in retaining meaningful relationships within complex datasets.
weaknesses: * The paper is about visualization. But there are no figures that visualize any real-world datasets in the main paper. There are some visualizations in the appendix. However, the results are even worse than t-SNE and UMAP. Especially when the colors are removed in an unsupervised setting, the MERCAT results are messy and show little useful global pattern.
* The paper addresses faithful visualization. But faithfulness is not defined.
* The method is expensive O(n^3) and has to rely on subsampling data points.

---

confidence: 2
contribution: 2
initial_rating: 6
presentation: 2
questions: 1. Why choose these parts, i.e., Smiley, Mammonth, and Circle in the embeddings of low-dimensional examples？The related reasons can be given here.

2. Why choose PCA instead of the others as the manner of reduction for robustness, whihc is shown in Algorithm 1.
soundness: 2
strengths: 1. Originality. he authors suggest a new perspective on LDE learning, reconstructing angles between data points. Then they show that this approach, MERCAT, yields good reconstruction across a diverse set of experiments and metrics, and preserve structures well across all scales. 

2. Quality and Clarity. There are no technical errors, and the presentation and writing are clear.

3. Significance. The authors show that this approach, MERCAT, yields good reconstruction across a diverse set of experiments and metrics, and preserve structures well across all scales. Compared to existing work, the proposed approach also has a simple formulation, facilitating future theoretical analysis and algorithmic improvements.
summary: The current generation of LDE approaches focus on reconstructing local distances between any pair of samples correctly, often outperforming traditional approaches aiming at all distances. For these approaches, global relationships are, however, usually strongly distorted, often argued to be an inherent trade-off between local and global structure learning for embeddings. The authors suggest a new perspective on LDE learning, reconstructing angles between data points. Then they show that this approach, MERCAT, yields good reconstruction across a diverse set of experiments and metrics, and preserve structures well across all scales. Compared to existing work, the proposed approach also has a simple formulation, facilitating future theoretical analysis and algorithmic improvements.
weaknesses: I am absolutely not in this field and the comments from me are not relatively professional. The comments in the following are just raised from the presentation or organization.

1. In the abstract part, the authors just state that their approach also has a simple formulation, facilitating future theoretical analysis and algorithmic improvements compared to existing work. I think the authors can add the specific values of performance gains here to better show the effectiveness of the proposed method.

2. The authors report the angle preservation, distance preservation, neighborhood preservation and density preservation between computed low-dimensional embeddings and original data. However, the authors use "-" to indicate method do not terminate within 24h. As far as I know, some experiments still report the final results that do not terminate within 24h. I think the authors can give more detailed explanations here.

3. For the real data results shown in Table 1, the authors should add more recent methods for comparison in the experiment to better show the effectiveness of the proposed method.

4. In the related work part, the authors just discuss three methods proposed after 2023, then the authors can add more recent related work for discussion in this part.

Meta-review: The paper proposes an algorithm for low-dimensional embedding (LDE) that aims to preserve angles between data points. Specifically, the high-dimensional data are mapped onto the 2-dimensional unit sphere. The proposed method is illustrated by numerical experiments.

Reviewers generally appreciate that the proposed idea is interesting along with its theoretical analysis. However, there are concerns both about the formulation and numerical results (see below)

Decision: Reject