:
AABI 2024Under review for the Workshop at the 6th Symposium on Advances in
Approximate Bayesian Inference, 2024

# Towards One Model for Classical Dimensionality Reduction: A Probabilistic
Perspective on UMAP and t-SNE

Aditya Ravuri
University of Cambridge
Neil D. Lawrence
University of Cambridge

###### Abstract

This paper shows that the dimensionality reduction methods, UMAP and t-SNE,
can be approximately recast as MAP inference methods corresponding to a
generalized Wishart-based model introduced in Ravuri et al. (2023). This
interpretation offers deeper theoretical insights into these algorithms, while
introducing tools with which similar dimensionality reduction methods can be
studied.

##  1 Introduction

In the realm of single-cell biology and various other domains with complex,
high-dimensional data, dimensionality reduction (DR) algorithms are essential
tools for uncovering the underlying structure of data. These algorithms, which
include very widely-used techniques like t-SNE (van der Maaten and Hinton,
2008) and UMAP (McInnes et al., 2020), are especially valuable for visualizing
data manifolds, enabling downstream processing, and the discovery of
insightful patterns. A deeper comprehension of these algorithms and their
theoretical underpinnings is crucial for advancing their applicability,
particularly when prior information is available and improving their
interpretability. Our work builds on (and aims to unify fully) the ProbDR
framework, which interprets classical DR methods through a probabilistic lens
to enable the communication of assumptions, integration of prior knowledge and
accounting for noise and confounders.

Ravuri et al. (2023) introduced ProbDR as a framework with two main
interpretations: UMAP and t-SNE corresponding to inference over an adjacency
matrix, and other classical DR algorithms that utilize eigendecomposition of a
PSD matrix as inference using a Wishart generative model on the
covariance/precision matrix.

In this work, we further simplify the framework, moving away from the
variational interpretation and propose that all algorithms with ProbDR
interpretations (and hence most classical DR methods) can be written as MAP
inference algorithms given the model,

$\mathbf{S}|\mathbf{X}\sim\mathcal{W}^{\\{-1\\}}(\mathbf{XX}^{T}+\epsilon\mathbf{H}K_{t}(\mathbf{X})\mathbf{H}+\gamma\mathbf{I},\nu),$
(1)

where $\mathbf{S}\in S^{+}_{n}$ is an estimate of a covariance matrix,
$\mathbf{X}\in\mathbb{R}^{n,q}$ corresponds to the set of low ($q$)
dimensional latent variables, $\mathbf{H}=\mathbf{I}-\mathbf{11}^{T}/n$ is a
centering matrix, and
$K_{t}(\mathbf{X}_{i},\mathbf{X}_{j})=(1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2})^{-1}$
is the Student-t kernel.

We aim to interpret these DR algorithms as MAP inference methods instead of
variational methods as studied in Ravuri et al. (2023); Van Assel et al.
(2022). This unifies t-SNE and UMAP-like algorithms with the other algorithms
and provides semantic interpretation to UMAP and t-SNE. Additionally, we hope
that the tools introduced in this paper will provide researchers with more
machinery to understand the behaviour of latent variable models.

##  2 Background

The ProbDR framework showed that many classical DR algorithms can be expressed
as inference algorithms corresponding to a probabilistic model. Algorithms
that set the embedding of a high-dimensional dataset
$\mathbf{Y}\in\mathbb{R}^{n,d}$ in terms of the eigenvectors of a positive-
semi-definite matrix were shown in Ravuri et al. (2023) to correspond to a
two-step inference process, where,

  1. 1.

one first estimates a covariance matrix $\mathbf{S}(\mathbf{Y})$ or a
precision matrix $\mathbf{\Gamma}(\mathbf{Y})$ (which a graph Laplacian
$\mathbf{L}$ can be an estimate of),

  2. 2.

then estimates the embedding via maximum a-posteriori given one of the two
following models,

$\displaystyle\mathbf{S}|\mathbf{X}\sim\mathcal{W}\left(\mathbf{X}\mathbf{X}^{T}+\sigma^{2}\mathbf{I}_{n},d\right)\text{
or,}$
$\displaystyle\mathbf{\Gamma}|\mathbf{X}\sim\mathcal{W}\left((\mathbf{X}\mathbf{X}^{T}+\beta\mathbf{I}_{n})^{-1},d\right).$

PCoA, for example, is recovered using the first of these formulations with
$\mathbf{S}\equiv\mathbf{Y}\mathbf{Y}^{T}$. Setting $\epsilon=0$ in Equation 1
recovers these results.

In the case of UMAP and t-SNE, ProbDR did not specify a generative model for
either the data or the covariance but only a generative model for the
adjacency matrices that describe the nearest neighbour graph.

Specifically, ProbDR showed that the generative models corresponding to UMAP
and t-SNE can be seen as models for adjacency matrices $\mathbf{A}^{\prime}$,

$\small
p(\mathbf{A}^{\prime}|\mathbf{X})=\begin{cases}\text{Categorical}(\text{vec}(\mathbf{A}^{\prime})|w^{t}_{ij}(\mathbf{X}))&\text{t-SNE}\\\
\prod_{i>j}^{n}\text{Bernoulli}(\mathbf{A}^{\prime}_{ij}|w^{U}_{ij}(\mathbf{X}_{i},\mathbf{X}_{j}))&\text{UMAP}\end{cases}$
(2)

where,

$\displaystyle
w^{t}_{ij}=\dfrac{(1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2})^{-1}}{\sum_{k\neq
l}(1+\|\mathbf{X}_{k}-\mathbf{X}_{l}\|^{2})^{-1}},\text{ and
}w^{U}_{ij}=\dfrac{1}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}.$

The inference can be done as maximum a-posteriori estimation for $\mathbf{X}$,
given the binary, empirical, nearest-neighbour adjacency matrix
$\mathbf{A}^{\prime}(\mathbf{Y})_{ij}=\mathcal{I}(j\in N(i))$, where $N(i)$
represents the set of nearest neighbours of data point $\mathbf{Y}_{i}$. UMAP
and t-SNE are typically interpreted in a variational way, however the
inference trivially becomes MAP inference when we use
$\mathbf{A}^{\prime}(\mathbf{Y})$ as the variational data-dependent
distribution (due to the Ravuri et al. (2023), Appendix B.7, Lemma 13). This
interpretation is also presented in Damrich et al. (2022).

This simplification is reasonable due to the findings of Damrich and Hamprecht
(2021), where it was found that the relatively complex calculation of the
variational probabilities in t-SNE and UMAP can be replaced with simply the
adjacency matrices without loss of performance. Our initial experiments also
closely aligned with these findings.

Crucially, however, Becht et al. (2019); Damrich et al. (2022) note that the
optimisation process is equally as important. As part of an extensive study on
the nature of the t-SNE and UMAP loss functions, Damrich et al. (2022) then
show how the stochastic optimisation of t-SNE and UMAP can be interpreted to
be contrastive estimation with the loss

$\mathcal{L}(\mathbf{X})\propto-\mathbb{E}_{ij\sim
p}\log\left(\frac{w_{ij}(\mathbf{X})}{w_{ij}(\mathbf{X})+1}\right)-m\mathbb{E}_{ij\sim\xi}\log\left(1-\frac{w_{ij}(\mathbf{X})}{w_{ij}(\mathbf{X})+1}\right),$
(3)

where $p$ represents a discrete distribution that is uniform across the
nearest neighbour pairs and zero everywhere else. Similarly, $\xi$ represents
a uniform distribution over non-neighbours. $m$ (set to $2n_{-}$) is a
multiplicative hyperparameter proportional to the number of contrastive
negatives that affects the strength of repulsion. In Damrich et al. (2022),
$w_{ij}(\mathbf{X})=\frac{1}{\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}$ recovers
UMAP. This bound is important, as we found that a naive optimisation of the
Bernoulli likelihood in Equation 2 leads to a poor embedding (although
Appendix B offers more commentary on this bound as a Bernoulli likelihood and
provides more evidence for our claims).

For this work, we aim to work with such a contrastive loss function and
interpret it as a likelihood, but over the latents $\mathbf{X}$. This is
because we were particularly inspired by Nakamura et al. (2023), who showed
that contrastive learning methods could be seen as variational algorithms
(hence suggesting a link between t-SNE and contrastive learning) and by
Gutmann and Hyvärinen (2010), which shows that contrastive losses are
estimators of negative log-likelihoods. Damrich et al. (2022) also greatly
simplify the UMAP optimisation process.

##  3 Discussion

This section argues that inference with the generative model in Equation 1
approximately recovers UMAP and t-SNE-like algorithms.

###  3.1 The proposed model

We first consider MAP inference for $\mathbf{X}$ given the model,

$\displaystyle\nu\mathbf{L}|\mathbf{X}$
$\displaystyle\sim\mathcal{W}((-\alpha\mathbf{H}\mathbf{D}^{\prime}\mathbf{H}+\gamma\mathbf{I})^{-1},\nu),$
(4) $\displaystyle\mathbf{X}$
$\displaystyle\sim\text{Uniform}(-\infty,\infty),$

where $\mathbf{D^{\prime}}$ is a squared distance matrix, with elements
$\mathbf{D}^{\prime}_{ij}=\log(1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2})$. Let
$\mathbf{M}=-\mathbf{H}\mathbf{D}^{\prime}\mathbf{H}$. $\mathbf{M}$ is PSD
(with some interesting properties relating to isometric Euclidean embeddings,
see Section A.2). The graph Laplacian is computed as
$\mathbf{L}=(\mathbf{D-A})/\bar{d}$, with $\bar{d}$ being the average
degree111So that its inverse (i.e. the covariance) has diagonal elements
around one.. The adjacency matrix $\mathbf{A}_{ij}=1\text{ if }p(ij)\text{ (of
\lx@cref{creftype~refnum}{eqn:cne-bound}) }>0.$

Then, the log-likelihood given the model in Equation 4 is as follows (we focus
on just the data-dependent term),

$\displaystyle\log p(\mathbf{X}|\mathbf{L})$
$\displaystyle=-0.5\nu\text{tr}(\mathbf{L}(\alpha\mathbf{M}+\gamma\mathbf{I}))+0.5\nu\log\det(\alpha\mathbf{M}+\gamma\mathbf{I})+c$
$\displaystyle\propto-\alpha\text{tr}(\mathbf{L}\mathbf{M})+\log\det(\alpha\mathbf{M}+\gamma\mathbf{I})+k$
$\displaystyle=\alpha\text{tr}(\mathbf{L}\mathbf{H}\mathbf{D^{\prime}}\mathbf{H})+\log\det(\alpha\mathbf{M}+\gamma\mathbf{I})+k$
$\displaystyle=\alpha\text{tr}\left(\mathbf{L}\mathbf{D}^{\prime}\right)+\log\det(\alpha\mathbf{M}+\gamma\mathbf{I})+k\qquad\text{(trace
cyclic and }\mathbf{L}\text{ centered)}$
$\displaystyle=-\alpha\text{tr}(\mathbf{A}\mathbf{D}^{\prime})+\log\det(\alpha\mathbf{M}+\gamma\mathbf{I})+k\qquad\text{tr}(\mathbf{D^{\prime}})=0$
$\displaystyle=-\alpha\sum_{i}\mathbf{A}_{i}^{T}\mathbf{D}^{\prime}_{i}+\log\det(\alpha\mathbf{M}+\gamma\mathbf{I})+k$
$\displaystyle=-\alpha\sum_{i}\sum_{j}a_{ij}\log(1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2})+\log\det(\alpha\mathbf{M}+\gamma\mathbf{I})+k$

$\displaystyle=\alpha\dfrac{n_{+}}{\bar{d}}\mathbb{E}_{ij\sim
p}\log\left(\dfrac{1}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j})\|^{2}}\right)+\log\det(\alpha\mathbf{M}+\gamma\mathbf{I})+k$
$\displaystyle=\alpha n\mathbb{E}_{ij\sim
p}\log\left(\dfrac{w_{ij}(\mathbf{X}_{i},\mathbf{X}_{j})}{1+w_{ij}(\mathbf{X}_{i},\mathbf{X}_{j})}\right)+\log\det(\alpha\mathbf{M}+\gamma\mathbf{I})+k,$

with $n_{+}=n\bar{d}\approx 1.5nn_{\\#}$ being the number of total number of
edges and $n_{\\#}$ is the number of neighbours per point (typically set to
15). An important note here is that the model is misspecified (for example,
the variance implied by the covariance parameter and the data are quite
different, with the variance of the Wishart being much lower than the data
estimate).

Therefore, by minimising the first term of Equation 3, we maximise the data-
dependent term of the likelihood of Equation 4, and so, Equation 4 defines a
model for dimensionality reduction that in some ways is similar to t-SNE and
UMAP-like algorithms. Figure 1 shows embeddings of 30,000 digits from the
MNIST dataset obtained using this model. We also run this model on a suite of
other datasets (from Pedregosa et al. (2011); Deng (2012); sta ; Krumsiek et
al. (2011)) to show that we recover roughly similar embeddings as minimisation
of Equation 3 in Figure 2.

Figure 1: Embeddings of 30k MNIST digits obtained using inference within
Equation 4. Left: A run with both $\alpha=1/50$ and $\gamma=1/5$ treated as
hyper-parameters. Right: A run with $\alpha=1/50$ treated as a hyperparameter
and $\gamma$ optimised using maximum likelihood inference. Appendix B sheds
light into why this choice of hyper-parameters is performant, leading to the
embedding of Figure 4.

Figure 2: Left: embeddings on various datasets obtained using optimisation of
the CNE bound (Equation 3) compared with right: inference results using our
model.

.

###  3.2 CNE bounds can sometimes be approximated using the proposed model

Consider the negative CNE loss,

$\displaystyle-\mathcal{L}$ $\displaystyle\propto\dfrac{\alpha\nu
n}{2}\mathbb{E}_{ij\sim
p}\log\left(\frac{w_{ij}(\mathbf{X})}{w_{ij}(\mathbf{X})+1}\right)+\dfrac{\alpha\nu
n}{2}m\mathbb{E}_{ij\sim\xi}\log\left(1-\frac{w_{ij}(\mathbf{X})}{w_{ij}(\mathbf{X})+1}\right)$
$\displaystyle=-\dfrac{\nu}{2}\text{tr}(\mathbf{L}(\alpha\mathbf{M}+\gamma\mathbf{I}))+\dfrac{\alpha\nu
n}{2}m\mathbb{E}_{ij\sim\xi}\log\left(\frac{\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}\right)$
Sec. 3.1

Let the first term be represented as $\mathcal{T}_{a}$. We first approximate
the distances $\mathbf{D}^{\prime}_{ij}=\log(1+d_{ij}^{2})\approx
0.5(1+d_{ij}^{2}-1/(1+d_{ij}^{2}))$ within $\mathcal{T}_{a}$ as this is a good
approximation for the log1p function near 0. Therefore,

$\displaystyle\mathbf{M}\approx-\mathbf{H}(\mathbf{11^{T}}+0.5\mathbf{D}^{2}-0.5K_{t})\mathbf{H}=\mathbf{X}\mathbf{X}^{T}+\mathbf{H}K_{t}\mathbf{H}.$
App. A.1

This approximation sheds light on the behaviour of the covariance parameter of
Equation 4 and draws the link to Equation 1.

Then,

$\displaystyle-\mathcal{L}$
$\displaystyle\propto\mathcal{T}_{a}+\dfrac{\alpha\nu
n}{2}m\mathbb{E}_{ij\sim\xi}\log\left(\frac{\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}\right)$
$\displaystyle\approx\mathcal{T}_{a}-\dfrac{\alpha\nu
mn}{2}\mathbb{E}_{ij\sim\xi}\left(\frac{1}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}\right)\qquad(\text{for
large distances})$ $\displaystyle\approx\mathcal{T}_{a}-\dfrac{\alpha\nu
mn}{2}\mathbb{E}_{ij\sim\mathcal{U}}\left(\frac{1}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}\right)\qquad\text{as
}n>>n_{\\#}$ $\displaystyle=\mathcal{T}_{a}-\dfrac{\alpha\nu
mn}{2n^{2}}\sum_{ij}\frac{1}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}$
$\displaystyle=\mathcal{T}_{a}+\dfrac{\alpha\nu
m}{2}\text{tr}\left(-\dfrac{1}{n}\mathbf{11^{T}}K_{t}-\dfrac{1}{n}K_{t}\mathbf{11^{T}}+\dfrac{1}{n^{2}}\mathbf{11^{T}}K_{t}\mathbf{11^{T}}\right)$
$\displaystyle=\mathcal{T}_{a}+\dfrac{\nu
m}{2}\text{tr}\left(\alpha\mathbf{H}K_{t}\mathbf{H}\right)+k$
$\displaystyle=-\dfrac{\nu}{2}\text{tr}(\mathbf{L}(\alpha\mathbf{M}_{\text{aprx}}+\gamma\mathbf{I}))+\dfrac{\nu}{2}\text{tr}\left(\alpha
m\mathbf{M}_{\text{aprx}}+\gamma m\mathbf{I}\right)-\dfrac{\alpha
m\nu}{2}\text{tr}(\mathbf{X}\mathbf{X}^{T})+k$
$\displaystyle\propto-\dfrac{m\nu}{2}\text{tr}(m^{-1}\mathbf{L}(\dfrac{\alpha}{\gamma}\mathbf{M}_{\text{aprx}}+\mathbf{I}))+\dfrac{m\nu}{2}\text{tr}\left(\dfrac{\alpha}{\gamma}\mathbf{M}_{\text{aprx}}+\mathbf{I}\right)-\dfrac{\alpha
m\nu}{2\gamma}\text{tr}(\mathbf{X}\mathbf{X}^{T})+k$
$\displaystyle\approx-\dfrac{m\nu}{2}\text{tr}(m^{-1}\mathbf{L}(\dfrac{\alpha}{\gamma}\mathbf{M}_{\text{aprx}}+\mathbf{I}))+\dfrac{m\nu}{2}\log\det\left(\dfrac{\alpha}{\gamma}\mathbf{M}_{\text{aprx}}+\mathbf{I}\right)-\dfrac{\alpha
m\nu}{2\gamma}\text{tr}(\mathbf{X}\mathbf{X}^{T})+k\quad(\text{large }\gamma)$
$\displaystyle=\log\mathcal{W}(\nu\mathbf{L}|m^{-1}(\alpha\mathbf{M}_{\text{aprx}}+\gamma\mathbf{I})^{-1},m\nu)+\log\mathcal{N}(\mathbf{X}|\mathbf{0},\gamma\mathbf{I}/\alpha
m\nu)$

This shows the NCE bound can be approximated by the likelihood of our proposed
model with a diffuse prior on $\mathbf{X}$. Note that the model is still
misspecified. A comparison of embeddings used with this approximation and the
CNE bound are shown in Figure 3.

Figure 3: Left: Embedding of 30k digits from MNIST obtained using the CNE
bound and right: using our approximation. For this, we set $n_{-}=1$, which
was the best performing case.

##  4 Conclusion

We present a probabilistic interpretation of t-SNE and UMAP-like algorithms,
showing that they correspond to inference within misspecified generative
models for the covariance/precision with fixed scale parameters, and a choice
of a covariance that describes non-linear functions within a Gaussian process
context. We hope that this serves as a foundation to further refine these
models (particularly based on their regularisation terms) such that they
emulate results from t-SNE and UMAP-like algorithms.

Figure 4: Embedding of 30k MNIST digits resulting from ProbDR inference using
a scaled Student-t kernel. See Appendix B.

AR would like to thank Francisco Vargas for helpful discussions and a
studentship from the Accelerate Programme for Scientific Discovery.

## References

  * (1) Statlog (Shuttle).  UCI Machine Learning Repository.  URL https://doi.org/10.24432/C5WS31.
  * Becht et al. (2019) Etienne Becht, Leland McInnes, John Healy, Charles-Antoine Dutertre, Immanuel W. H. Kwok, Lai Guan Ng, Florent Ginhoux, and Evan W. Newell.  Dimensionality reduction for visualizing single-cell data using UMAP.  _Nature Biotechnology_ , 37(1):38–44, Jan 2019.  ISSN 1546-1696.  10.1038/nbt.4314.  URL https://doi.org/10.1038/nbt.4314.
  * Bhatia (2007) Rajendra Bhatia.  _Positive Definite Matrices_.  Princeton University Press, 2007.  ISBN 9780691129181.  URL http://www.jstor.org/stable/j.ctt7rxv2.
  * Damrich and Hamprecht (2021) Sebastian Damrich and Fred A Hamprecht.  On umap's true loss function.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, _Advances in Neural Information Processing Systems_ , volume 34, pages 5798–5809. Curran Associates, Inc., 2021.  URL https://proceedings.neurips.cc/paper_files/paper/2021/file/2de5d16682c3c35007e4e92982f1a2ba-Paper.pdf.
  * Damrich et al. (2022) Sebastian Damrich, Niklas Böhm, Fred A Hamprecht, and Dmitry Kobak.  From $t$-sne to umap with contrastive learning.  In _The Eleventh International Conference on Learning Representations_ , 2022.
  * Deng (2012) Li Deng.  The mnist database of handwritten digit images for machine learning research [best of the web].  _IEEE Signal Processing Magazine_ , 29(6):141–142, 2012.  10.1109/MSP.2012.2211477.
  * Gupta and Nagar (2018) Arjun K Gupta and Daya K Nagar.  _Matrix variate distributions_.  Chapman and Hall/CRC, 2018.
  * Gutmann and Hyvärinen (2010) Michael Gutmann and Aapo Hyvärinen.  Noise-contrastive estimation: A new estimation principle for unnormalized statistical models.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pages 297–304. JMLR Workshop and Conference Proceedings, 2010.
  * Khare (2019) Apoorva Khare.  Schoenberg: from metric geometry to matrix positivity, Apr 2019.  URL https://math.iisc.ac.in/seminar-slides/2019/2019-04-12-ApoorvaKhare.pdf.
  * Krumsiek et al. (2011) Jan Krumsiek, Carsten Marr, Timm Schroeder, and Fabian J Theis.  Hierarchical differentiation of myeloid progenitors is encoded in the transcription factor network.  _PLoS One_ , 6(8):e22649, August 2011.
  * McInnes et al. (2020) Leland McInnes, John Healy, and James Melville.  UMAP: Uniform manifold approximation and projection for dimension reduction, 2020.
  * Nakamura et al. (2023) Hiroki Nakamura, Masashi Okada, and Tadahiro Taniguchi.  Representation uncertainty in self-supervised learning as variational inference.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 16484–16493, 2023.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011.
  * Ravuri et al. (2023) Aditya Ravuri, Francisco Vargas, Vidhi Lalchand, and Neil D Lawrence.  Dimensionality reduction as probabilistic inference.  In _Fifth Symposium on Advances in Approximate Bayesian Inference_ , 2023.  URL https://arxiv.org/pdf/2304.07658.pdf.
  * Schoenberg (1935) I. J. Schoenberg.  Remarks to maurice frechet’s article “sur la definition axiomatique d’une classe d’espace distances vectoriellement applicable sur l’espace de hilbert.  _Annals of Mathematics_ , 36(3):724–732, 1935.  ISSN 0003486X.  URL http://www.jstor.org/stable/1968654.
  * Van Assel et al. (2022) Hugues Van Assel, Thibault Espinasse, Julien Chiquet, and Franck Picard.  A probabilistic graph coupling view of dimension reduction.  _Advances in Neural Information Processing Systems_ , 35:10696–10708, 2022.
  * van der Maaten and Hinton (2008) Laurens van der Maaten and Geoffrey Hinton.  Visualizing data using t-SNE.  _Journal of Machine Learning Research_ , 9(86):2579–2605, 2008.  URL http://jmlr.org/papers/v9/vandermaaten08a.html.

##  Appendix A Centered Distance Based Results

###  A.1 The MDS Matrix is a Gram matrix

Assume centered $\mathbf{X}$ (i.e. the column means of $\mathbf{X}$ are zero).
Then,

$\displaystyle-\mathbf{H}\mathbf{D}\mathbf{H}$
$\displaystyle=-\left(\mathbf{I}-\frac{1}{n}\mathbf{11}^{T}\right)\mathbf{D}\left(I-\frac{1}{n}\mathbf{11}^{T}\right)$
$\displaystyle=-\left(\mathbf{I}-\frac{1}{n}\mathbf{11}^{T}\right)\left(\tilde{\mathbf{D}}\mathbf{11}^{T}+\mathbf{11}^{T}\tilde{\mathbf{D}}-2\mathbf{X}\mathbf{X}^{T}\right)\left(I-\frac{1}{n}\mathbf{11}^{T}\right)$
$\displaystyle...$
$\displaystyle=2\mathbf{X}\mathbf{X}^{T}-\frac{2}{n}\mathbf{X}\mathbf{X}^{T}\mathbf{11}^{T}-\frac{2}{n}\mathbf{11}^{T}\mathbf{X}\mathbf{X}^{T}+\frac{2}{n^{2}}2\mathbf{11}^{T}\mathbf{X}\mathbf{X}^{T}\mathbf{11}^{T}$
$\displaystyle=2\mathbf{X}\mathbf{X}^{T}$

###  A.2 -$\mathbf{H}\mathbf{D}^{\prime}\mathbf{H}$ is PSD

$\mathbf{D}_{ij}^{{}^{\prime}0.5}$ (and $\mathbf{D}^{\prime}_{ij}$) are valid
distance metrics. We only need to show that $-\mathbf{D^{\prime}}$ is CPSD,
then, there exists an isometric Euclidean embedding Khare (2019); Schoenberg
(1935).
$-\mathbf{D^{\prime}}=\log\dfrac{1}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}$.
Note that the inner term is a kernel, hence is PSD. The $\log$ function
preserves CPSD-ness (Bhatia, 2007).

##  Appendix B A different interpretation of Equation 3

If the negative samples coefficient of the second term of Equation 3 is
absorbed into the logarithm, and the first term ignored, this implies that the
edge probability is,

$\displaystyle\mathbb{P}(\mathbf{A}_{ij}=1)$
$\displaystyle=1-\left(1-\dfrac{1}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}\right)^{1.5n_{\\#}n_{-}/n}$
$\displaystyle\approx\dfrac{1.5n_{\\#}n_{-}}{n}\log\left(1+\dfrac{1}{\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}+\epsilon}\right).$

We further approximate this probability as the following based on empirical
observations that it did not seem to affect the embeddings negatively,

$\displaystyle\mathbb{P}(\mathbf{A}_{ij}=1)$
$\displaystyle\approx\dfrac{1.5n_{\\#}n_{-}/n}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}.$

Using this definition, we estimate the parameter of a Wishart distribution
over the normalized graph Laplacian $\mathbf{L}$ by moment matching
(specifically such that the functional terms of the Wishart variance match the
Bernoulli distribution’s variances). This leads to the covariance estimate,

$\mathbf{M}_{ij}=\dfrac{\mathbf{p}_{ij}}{\sum_{a}\mathbf{p}_{ia}\sum_{b}\mathbf{p}_{ib}},$

where
$\mathbf{p}_{ij}=\dfrac{1.5n_{\\#}n_{-}/n}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}$.
It’s easy to see that this is PD. Plugging this in place of $\alpha\mathbf{M}$
into our model in Equation 4 leads to the embedding in LABEL:fig:omp-est. Note
that the coefficient is around $0.01$ when $n=10k$, which might shed light
into the choice for $\alpha$ that is most performant in Figure 4.

##  Appendix C A Matrix-t Perspective

The density of the matrix-t distribution, with the column covariance set to
$\mathbf{I}$ is expressed as (Gupta and Nagar, 2018),

$\displaystyle\log p(\mathbf{Y}|\Sigma)$
$\displaystyle=-(\alpha+n/2)\log\left|\mathbf{I}_{n}+\dfrac{\beta}{2}\Sigma^{-1}\mathbf{Y}\mathbf{Y}^{T}\right|-\dfrac{d}{2}\log|\Sigma|+c.$

Let $\mathbf{S}=\mathbf{\Gamma}^{-1}$ be an invertible estimator of
$\mathbf{Y}\mathbf{Y}^{T}$. Note that,

$\displaystyle\log\left|\mathbf{I}_{n}+\dfrac{\beta}{2}\Sigma^{-1}\mathbf{S}\right|$
$\displaystyle=\log\left|\dfrac{\beta}{2}\Sigma^{-1}\mathbf{S}\right|\left|\dfrac{2}{\beta}\mathbf{\Gamma}\Sigma+\mathbf{I}_{n}\right|$
$\displaystyle=\log\left|\dfrac{2}{\beta}\mathbf{\Gamma}\Sigma+\mathbf{I}_{n}\right|-\log\left|\Sigma\right|+k.$

Therefore the density can be rewritten in terms of the precision matrix
$\mathbf{\Gamma}$,

$\displaystyle\log p(\mathbf{Y}|\Sigma)$
$\displaystyle=\dfrac{2\alpha+n-d}{2}\log|\Sigma|-(\alpha+n/2)\log\left|\dfrac{2}{\beta}\mathbf{\Gamma}\Sigma+\mathbf{I}_{n}\right|+c.$

For reference, the negative log-density of a matrix Cauchy distribution can be
written as follows (adapted from ),

$\displaystyle\mathcal{L}_{t_{1}}(\Sigma)=\dfrac{d+n}{2}\log|I+\mathbf{L}\Sigma|-\dfrac{n}{2}\log|\Sigma|+c.$

The first term of the CNE bound can also be approximated as,

$\displaystyle-\mathbb{E}_{ij\sim
p}\log\left(\frac{w_{ij}(\mathbf{X})}{w_{ij}(\mathbf{X})+1}\right)$
$\displaystyle=-\mathbb{E}_{ij\sim
p}\log\left(\frac{1}{1+\frac{1}{w_{ij}(\mathbf{X})}}\right)$
$\displaystyle=-\mathbb{E}_{ij\sim
p}\log\left(\frac{1}{2+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}\right)$
$\displaystyle=\mathbb{E}_{ij\sim
p}\log\left(2+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}\right)$
$\displaystyle\approx
n^{-2}\sum_{ij}a_{ij}\log\left(1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}\right)$
$\displaystyle\leq\log(2+n^{-2}\sum_{ij}a_{ij}\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2})$
$\displaystyle=\log(1+2n^{-2}\text{tr}(\mathbf{L}\mathbf{X}\mathbf{X}^{T}+\mathbf{I}))$
$\displaystyle\approx\log|\mathbf{I}+\mathbf{L}(\mathbf{X}\mathbf{X}^{T}+\mathbf{I})|.$

Following the methodology of Section 3.1 however leads to a very Laplacian-
Eigenmaps-like solution. The regularisation term is paramount to the quality
of embeddings obtained using such methods.

##  Appendix D An Analysis of the Adjacency Probabilities

The generative probabilities of t-SNE, in some sense, model the probability
with which an edge is the shortest edge on the entire distance graph,
representing the shortest edge. The probability of this is proportional to
$w_{ij}(\mathbf{X})=\frac{1}{1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}}$
(implied by the contrastive loss used in Damrich et al. (2022); this is also
evident from the construction of the algorithm).

In this section, we assume that there’s a latent dataset described by a
generative model $\mathbf{X}\rightarrow\mathbf{Y}$ that’s unobserved, but
whose statistics then affect the data adjacency matrix $\mathbf{A}$ through
$\mathbf{Y}\rightarrow\mathbf{A}$, that is observed.

Assuming that a generative model for data $\mathbf{Y}$ given latents
$\mathbf{X}$ exists, we show by simulation that
$\mathbb{P}(\mathcal{A}_{ij}(\mathbf{Y})=1|\mathbf{X})=\mathbb{P}\left(\mathcal{I}\left[\|\mathbf{Y}_{i}-\mathbf{Y}_{j}\|^{2}=\underset{k>l}{\text{argmin}}\|\mathbf{Y}_{k}-\mathbf{Y}_{l}\|^{2}\right]=1\mathbf{X}\right)$
cannot be proportional to the form assumed in t-SNE if a Gaussian process
model,

$\mathbf{Y}|\mathbf{X}\sim\mathcal{MN}(\mathbf{0},\mathbf{X}\mathbf{X}^{T}+\sigma^{2}\mathbf{I},\mathbf{I}),$

is assumed, but can be achieved by both matrix Cauchy distributiones with a
linear kernel and Gaussian processes with a kernel that is the sum of linear
and smooth kernels. This is intuitive, as for these probabilities to decay to
zero, in the Gaussian case, the kernel must be non-stationary. In the case of
a matrix Cauchy distribution, this is due to its extreme-value properties.

Monte-Carlo simulations were done where a matrix Cauchy distribution (left on
Figure 5) and a Gaussian process (right) with dot product kernels are used to
generate high-dimensional data, and hence distance matrices between the high-
dimensional data points. Then, the log probability that the $ij^{th}$ element
of the distance matrix is the smallest is computed and plotted against the
Euclidean distances of the corresponding $\mathbf{X}$ points. These show that
using a Gaussian process with a linear kernel produces probabilities that are
linear functions of the latent Euclidean distance. A matrix Cauchy prior on
the other hand, induces proximity probabilities of the right shape.

Figure 5: Monte-Carlo simulations showing the probability with which a
distance $d_{ij}(\mathbf{Y})$ is the minimum throughout the entire distance
matrix, plotted as a function of the Euclidean distance between the
$\mathbf{X}$.

As the linear kernel induces extremely small probabilities for larger
distances in $\mathbf{X}$, we also use importance sampling to simulate tail
probabilities and ensure that these behave linearly as a function of the
distances in $\mathbf{X}$. For this, we need a full joint distribution of
distances, which is approximated below.

######  Theorem D.1 (Distribution of normal distances).

Assume that $\mathbf{Y}$ is distributed as,

$\begin{bmatrix}\mathbf{Y}_{i}\\\
\mathbf{Y}_{j}\end{bmatrix}\sim\mathcal{MN}\left(\mu,\begin{bmatrix}k_{ii}&k_{ij}\\\
k_{ji}&k_{jj}\end{bmatrix},\mathbf{I}_{d}\right).$

Then, the following hold. Firstly, denoting
$d_{ij}^{2}=\|\mathbf{Y}_{i}-\mathbf{Y}_{j}\|^{2}$, the marginal distribution
is given by,

$d_{ij}^{2}\sim\Gamma\left(k=\dfrac{d}{2},\theta=2(k_{ii}+k_{jj}-2k_{ij})\right).$

As a consequence, $\mathbb{E}(d_{ij}^{2})=d*\tilde{k}_{ij}\text{ and
}\mathbb{V}(d_{ij}^{2})=2d*\tilde{k}_{ij}^{2}$, where
$\tilde{k}_{ij}=k_{ii}+k_{jj}-2k_{ij}$.

Additionally,

$\mathbb{C}(d_{ij}^{2},d_{mn}^{2})=2d*(k_{im}+k_{jn}-k_{in}-k_{jm})^{2}.$

This is a useful fact as the upper triangle of the distance matrix is
approximately normal due to the central limit theorem with increasing $d$.
Proved in Appendix E.

Results of this experiment are given in Figure 6, which confirms that a
Gaussian process with a linear (dot product) kernel produces probabilities
that are approximately linear functions of the latent distance, which are
unlike the assumed tSNE adjacency probabilities.

Figure 6: Illustration of $\mathbb{P}(argmin(d_{ij}^{2})=ij)$ through
importance sampling for a linear GP.

Along different lines, we can think of the UMAP generative model as a model
for adjacencies
$\mathbf{A}_{ij}\sim\text{Bernoulli}(1/(1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}))$
as stated in ProbDR. Using just the mean the variance in Theorem D.1, we can
find the parameters of a covariance
$\mathbf{S}=\alpha\mathbf{X}\mathbf{X}^{T}+\beta
K_{t}(\mathbf{X})+\gamma\mathbf{I}$, using a black-box optimiser such that the
probability of a pair of points $ij$ having a distance smaller than
$\epsilon$,

$\mathbb{P}(d_{ij}^{2}<\epsilon)\approx\Phi\left(\dfrac{\epsilon-d\tilde{k}_{ij}}{\sqrt{2d}\tilde{k}_{ij}}\right),$

with
$\tilde{k}_{ij}=\alpha\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}+2\gamma+2\beta(1-1/(1+\|\mathbf{X}_{i}-\mathbf{X}_{j}\|^{2}))$,
is as in UMAP. This results in the parameters $\alpha\approx 0.1,\beta\approx
0.3,\gamma\approx 0.6$, with the fit being near-perfect. However, we reason
that this is quite a coarse approximation and a naive application of this
covariance within ProbDR results in a mediocre embedding.

##  Appendix E An Approximate Joint Distribution of Euclidean Distances

The first part of the theorem is given in Ravuri et al. 2023, reproduced
below.

$\displaystyle\forall k:d^{\prime}_{ij}\equiv y_{i}^{k}-y_{j}^{k}$
$\displaystyle\sim\mathcal{N}(0,k_{ii}+k_{jj}-2k_{ij})\overset{d}{=}\sqrt{k_{ii}+k_{jj}-2k_{ij}}Z$
$\displaystyle\Rightarrow
d_{ij}^{2}\equiv\|\mathbf{Y}_{i}-\mathbf{Y}_{j}\|^{2}=\sum_{k}^{d}(y^{k}_{i}-y^{k}_{j})^{2}$
$\displaystyle\overset{d}{=}(k_{ii}+k_{jj}-2k_{ij})\sum_{k}^{d}Z_{k}^{2}$
$\displaystyle\overset{d}{=}(k_{ii}+k_{jj}-2k_{ij})\chi^{2}_{d}$
$\displaystyle\overset{d}{=}\Gamma(k=d/2,\theta=2(k_{ii}+k_{jj}-2k_{ij})).$

The covariance between $d_{ij}^{2}$ and $d_{mn}^{2}$ can be computed as
follows. Let,

$\displaystyle
d^{\prime}_{ij}=\mathbf{Y}_{id}-\mathbf{Y}_{jd}\;\;\text{and}\;\;d^{\prime}_{mn}=\mathbf{Y}_{md}-\mathbf{Y}_{nd}.$

We can then derive some important moments as follows,

$\mathbb{E}(d^{\prime}_{ij})=0,$
$\mathbb{V}(d^{\prime}_{ij})=\mathbb{E}(d^{{}^{\prime}2}_{ij})=k_{ii}+k_{jj}-2k_{ij},$

$\displaystyle\mathbb{C}(d^{\prime}_{ij},d^{\prime}_{mn})$
$\displaystyle=\mathbb{C}(Y_{id}-Y_{jd},Y_{md}-Y_{nd})$
$\displaystyle=\mathbb{C}(Y_{id},Y_{md})-\mathbb{C}(Y_{id},Y_{nd})-\mathbb{C}(Y_{jd},Y_{md})+\mathbb{C}(Y_{jd},Y_{nd})$
$\displaystyle=k_{im}+k_{jn}-k_{in}-k_{jm},$

Then,

$\displaystyle\mathbb{C}(d_{ij}^{2},d_{mn}^{2})$
$\displaystyle=\mathbb{C}\left(\sum_{d_{1}}(\mathbf{Y}_{id_{1}}-\mathbf{Y}_{jd_{1}})^{2},\sum_{d_{2}}(\mathbf{Y}_{md_{2}}-\mathbf{Y}_{nd_{2}})^{2}\right)$
$\displaystyle=\sum_{d_{1}}\sum_{d_{2}}\mathbb{C}((\mathbf{Y}_{id_{1}}-\mathbf{Y}_{jd_{1}})^{2},(\mathbf{Y}_{md_{2}}-\mathbf{Y}_{nd_{2}})^{2})$
linearity
$\displaystyle=\sum_{d}\mathbb{C}(d^{{}^{\prime}2}_{ij},d^{{}^{\prime}2}_{mn})$
independence
$\displaystyle=d*\mathbb{C}(d^{{}^{\prime}2}_{ij},d^{{}^{\prime}2}_{mn})$
$\displaystyle=d*\left[\mathbb{E}[d^{{}^{\prime}2}_{ij}d^{{}^{\prime}2}_{mn}]-\mathbb{E}[d^{{}^{\prime}2}_{ij}]\mathbb{E}[d^{{}^{\prime}2}_{mn}]\right]$
$\displaystyle=d*\left[\mathbb{E}(d^{{}^{\prime}2}_{ij})\mathbb{E}(d^{{}^{\prime}2}_{mn})+2\mathbb{E}^{2}(d^{\prime}_{ij}d^{\prime}_{mn})-\mathbb{E}[d^{{}^{\prime}2}_{ij}]\mathbb{E}[d^{{}^{\prime}2}_{mn}]\right]$
Isserlis’ theorem
$\displaystyle=2d*\mathbb{E}^{2}(d^{\prime}_{ij}d^{\prime}_{mn})$
$\displaystyle=2d*(k_{im}+k_{jn}-k_{in}-k_{jm})^{2}.$