A Heat Diffusion Perspective on Geodesic Preserving Dimensionality Reduction
Authors: Guillaume Huguet
Alexander Tong
Edward De Brouwer
Yanlei Zhang
Guy Wolf
Ian Adelstein
Smita Krishnaswamy
Venue: NeurIPS 2023 poster (2023)
Keywords: manifold learning
heat diffusion
geodesic
metric preserving
dimensionality reduction
embedding
TL;DR: We propose a heat-geodesic dissimilarity for dimensionality reduction, and relate it to other methods.

Abstract: Diffusion-based manifold learning methods have proven useful in representation learning and dimensionality reduction of modern high dimensional, high throughput, noisy datasets. Such datasets are especially present in fields like biology and physics. While it is thought that these methods preserve underlying manifold structure of data by learning a proxy for geodesic distances, no specific theoretical links have been established. Here, we establish such a link via results in Riemannian geometry explicitly connecting heat diffusion to manifold distances. In this process, we also formulate a more general heat kernel based manifold embedding method that we call heat geodesic embeddings. This novel perspective makes clearer the choices available in manifold learning and denoising. Results show that our method outperforms existing state of the art in preserving ground truth manifold distances,  and preserving cluster structure in toy datasets. We also showcase our method on single cell RNA-sequencing datasets with both continuum and cluster structure, where our method enables interpolation of withheld timepoints of data. Finally, we show that parameters of our more general method can be configured to give results similar to PHATE (a state-of-the-art diffusion based manifold learning method) as well as SNE (an attraction/repulsion neighborhood based method that forms the basis of t-SNE).

A Heat Diffusion Perspective on Geodesic Preserving Dimensionality Reduction
Guillaume Huguet!*
Alexander Tong
Edward De Brouwer
Yanlei Zhang'
Guy Wolf
Ian Adelstein
Smita Krishnaswamy
Université de Montréal; Mila
Quebec Al Institute
Yale University
Abstract
Diffusion-based manifold learning methods have proven useful in representation learning and dimensionality reduction of modern high dimensional, high through- put, noisy datasets. Such datasets are especially present in fields like biology and physics. While it is thought that these methods preserve underlying manifold struc- ture of data by learning a proxy for geodesic distances, no specific theoretical links have been established. Here, we establish such link via results in Riemannian geometry explicitly connecting heat diffusion to manifold distances. this process, we also formulate more general heat kernel based manifold embedding method that we call heat geodesic embeddings. This novel perspective makes clearer the choices available in manifold learning and denoising. Results show that our method outperforms existing state-of-the-ard in preserving ground truth manifold distances, and preserving cluster structure in toy datasets. We also showcase method single cell ReA-sequeniing datasets with both continuum and cluster structure, where method enables interpolation of withheld timepoints of data. Finally, we show that parameters of our more general method can be configured to give results similar to PHATE (a state-of-the-art diffusion based manifold learning method) well as SNE (an attraction/repulsion neighborhood based method that forms the basis of t-SNE).
Introduction
The advent of high throughput and high dimensional data in various fields of science have made dimensionality reduction and visualization techniques an indispensable part of exploratory analysis. Diffusion-basee manifold learning methods, based the data diffusion operator, first defined in 5. have proven especially useful due to their ability to handle noise and density variations while preserving structure. As result, diffusion-basee
dimensionality reduction methods, such as PHATE (22. T-PHATE (3), and diffusion maps (5). have emerged as methods for analyzing high throughput noisy data in various situations. While these methods surmised to learn manifold geodesic distances, no specific theoretical links have been established. Here, we establish such a link by using Varadhan's formula (34 and a parabolic Harnack inequality [17 (2). which relate manifold distances to heat diffusion directly. This lens gives new insight into existing dimensionality reduction methods, including when they preserve geodesics, and suggests a new method for dimensionality reduction to explicitly preserve geodesics, which we call heat geodesic embedding Furthermore, based on our understanding of other methods 22 we introduce theoretically justified parameter choices that
"Equal contribution
I Co--senior authors
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

allow our method to have greater versatility in terms of distance denoising and emphasis on local versus global distances.
Generally, data diffusion operators are created by first computing distances between datapoints, transforming these distances into affinities by pointwise application of kernel function (like Gaussian kernel), and then row normalizing with or without first applying degree normalization into a Markovian diffusion operator P [ 14 33 The entries of P(x,)) y) then contain probabilities of diffusing (or random walk probabilities) from one datapoint to another. Diffusion maps and PHATE use divergences between these diffusion or random walk-based probability distributions P(x,-) and -(y.-) -) to design a diffusion-basee framework directly utilizes a heat kernel based distance, and offers a more comprehensive perspective to study these diffusion methods. By configuring parameters in our framework, we show how we can navigate a continuum of embeddings methods from PHATE 22 to Stochastic Neighbor Embedding (SNE) 11
distance that may not directly relate to manifold distance. Our
In summary, our contributions are as follows:
We define the heat-geodesic dissimilarity based on Varadhan's formula and the two-sided heat kernel bounds.
Based on this dissimilarity, we present a versatile geodesic--resserinng method for dimensionality reduction which we call heat geodesic embedding.
We establish a relationship between diffusion-based distances and the heat-geodesic dissimilarity.
We establish connections between our method and popular dimensionality reduction techniques such PHATE and SNE, shedding light their geodesic preservation and denoising properties based on modifications of the computed dissimilarity and distance preservation losses.
We empirically demonstrate the advantages of Heat Geodesic Embedding in preserving manifold geodesic distances in several experiments showcasing more faithful manifold distances in the embedding space, as well as our ability to interpolate data within the manifold.
Swiss Roll
UMAP
E-SNE
PHATE
Diffusion Maps
HeatGeo
Surches
Figure I: Embeddings of the Swiss roll (top) and Tree (bottom) datasets for different manifold learning methods. Our HeatGeo method correctly unrolls the Swiss roll while t-SNE and UMAP create undesirable artificial clusters.
Preliminaries
First, we introduce fundamental notions that form the basis of our manifold learning methods: Varad- han's formula R4] on manifold, diffusion processes graphs, efficient heat kernel approximations, and multidimensional scaling 4 12 16)
Varadhan's 's formula Varadhan's 's formula is a powerful tool in differential geometry that establishes
a connection between the heat kernel and the shortest path (geodesic) distance on a Riemannian manifold. Its versatility has led to widespread applications
machine learning 6 I0 [5] 27
29 Let (M,g) be closed Riemannian manifold, and A the Laplace-Beltrant operator on M. The heat kernel (a, y) on M is the minimal positive fundamental solution of the heat equation de Au with initial condition ho(2, y) (y). In a d-dimensional Euclidean space the heat -4/2 kernel is (x, y) (4Tt) 4/2 dd((,,)) /4c so that 44 log he( y) = 2dt log(4nt) + df(2, y) and

---

we observe the following limiting behavior:
lim 4t log (e = (0, y). 200
Varadhan B4 (see also (0 proved that eq. [ (now Varadhan's 's formula) holds more generally on complete Riemannian manifolds M, where d(x.y) is the geodesic distance on M, and the convergence is uniform over compact subsets of M. A related result for complete Riemannian manifolds that satisfy the parabolic Harnack inequality (which includes convex domains in Euclidean space and Riemannian manifolds with non-negative Ricci curvature) is the two-sided heat kernel bound 25 17], showing that for any e E 10 there exist constants c(t) and C(e) such that
c(E) exp V(x, vt)
(- d(r, y) 4((++e)t)
C(e) S h, & V(a,, Vt)
exp (- d(z, y) 4(1 c)(t (<((1-e))
(2)
where V(x, vt) is the volume of ball or radius vi centered at We denote this relation by (a y) = V(x, vf) exp(-d(<., y) /t) and note that it again recovers eq. in the t 0 limit, which is unsurprising as Varadhan's result holds more generally. More important for our purposes is that (x(2,)) ~ V(z, vt)- exp(-d(x; /)/ holds for > which will allow us to approximate geodesic distances d(zz y) from diffusion based estimation of the heat kernel hy he y) and the volume (zz, vt). In appendix C.3
we provide examples using inequality
Graph construction and diffusion Our construction starts by creating a graph from a point cloud dataset X of size n. We use a kernel function F x Rd + R such that the (weighted) adjacency matrix is Wij ;= x(x,, xÿ) for all o4, Iy E X. The kernel function could be a Gaussian kernel, constructed from a nearest neighbor graph. The resulting graph G is characterized by the set of nodes (an ordering of the observations), the adjaceney matrix, and the set of edges, i.e. pairs of nodes with non-zero weights. The graph Laplacian is an operator acting on signals on such that it mimics the negative of the Laplace operator. The combinatorial graph Laplacian matrix is defined as L == Q - W and its normalized version as L = In - Q--QQ/--00 1/1/2 1/2, where Q is a diagonal degree matrix with Qu := E3 W The Laplacian is symmetric positive semi-definite, and has an sigen-decomposition L = V' Throughout the presentation, we assume that Qa > for all i E {1, (l....,n,, nj. The Laplacian allows to define the heat equation on G, with respect to an initial signal fo E R on G:
i at +f = s.t. f(0) = fo EER
The solution of the above differential equation is obtained with the matrix exponential f(i) = e-L fo- and we define the heat kernel on the graph as H, == e-LL_ By eigendecompositton, we have H; -te-tApp TT The matrix H, is a diffusion matrix that characterizes how a signal propagate through the graph according to the heat equations.
Other diffusion matrices on graphs have also been investigated in the literature. The transition matrix P Q- W characterizing a random walk on the graph is another common diffusion matrix used for manifold learning such as PHATE 22 and diffusion maps [.] It is a stochastic matrix that converges to a stationary distribution π. Qu/ Σ, Qsi under mild assumptions.
Fast computation of heat diffusion Exact computation of the (discrete) heat kernel H, is com- putationally costly, requiring a full eigendecomposition in O(n²) time. Fortunately, multiple fast approximations have been proposed, including using orthogonal polynomials or the Euler backward methods. In this work, we use Chebyshev polynomials, as they have been shown to converge faster than other polynomials on this problem [[3].
Chebyshev polynomials are defined by the recursive relation [Te)KEN with Ta(y) 0, (v) and (1) 2y/Ta- (v) Tx-2(y) for A: > 2. Assuming that the largest eigenvalue is less than two (which holds for the normalized Laplacian), we approximate the heat kernel with the truncated polynomials of order K
K be Tx (L -I,)
== br,0
H. = PK (L,t))
2
k=
3

---

where the K + 1 scalar coefficients {br,i} depend on time and are evaluated with the Bessel function. Computing Px (L, requires K matrix-vector product and K + Bessel function evaluation. The expensive part of the computation are the matrix- products, which can be efficient if the Laplacian matrix is sparse. Interestingly, we note that the evaluation of T7 do not depend on the diffusion time. Thus, to compute multiple approximations of the heat kernel {px (L, t)) ()),))) LET. only necessitates reweighting the truncated polynomial {Tx kK)........ with the corresponding [7] sets of Bessel coefficients. The overall complexity is dominated by the truncated polynomial computation which takes O(K (E + n)) time where E' is the number non-zero values in L.
Another possible approximation is using the Euler backward method. It requires solving K systems of linear equations ft) (I, + using the Cholesky decomposition [O 28 We quantify the differences between the heat kernel
f0), which can be efficient for sparse matrices
approximations Appendix C
Metric multidimennsonal scaling Given dissimilarity function d between data points, metric multidimensional scaling (MDS) [16] finds an embedding & such that the difference between the given dissimilarity and the Euclidean distance in the embedded space is minimal across all data points. Formally, for a given function Rd x Rd + RT, MDS minimizes the following objective:
L(X) = Σ (d(a,, x.)) [,(x;;) ((x)))12) 1/2
(5)
In metric MDS the solution is usually found by the SMACOF algorithm [0), or stochastic gradient descent 37
3 Related Work
We review state-of-the-an embedding methods and contextualize them with respect to Heat Geodesic Embedding. A formal theoretical comparison of all methods is given Section 5 Given a set of high-dimensional datapoints, the objective of embedding methods is to create map that embeds the observations in lower dimensional space, while preserving distances or similarities. Different methods vary by their choice of distance dissimilarity functions, shown below.
Diffusion maps In diffusion maps [5), an embedding in & dimensions is defined via the first k non-trivial right eigenvectors of the t-steps random walk P weighted by their eigenvalues. The where ö;
embedding preserves the diffusion distance DMp(2,, æ,) == (, - ô, (1) vector such that (ō.), if j i and 0 otherwise, and π is the stationary distribution of P. Intuitively, DMp(2; XT) considers all the t-steps paths between 2; and zj. A larger diffusion time can seen as a low frequency graph filter, i.e. keeping only information from the low frequency transitions such has the stationary distributions. For this reason, using diffusion with t > helps denoising the relationship between observations.
This diffusion- based method preserves the potential distance (22 PHp |1- log 8, log 1/12. and justifies this approach using the log transformation to prevent nearest neighbors from dominating the distances. An alternative approach is suggested using a square root transformation. Part of our contributions is to justify the log transformation from a geometric point of view. The embedding is defined using multidimensional scaling, which we present below.
PHATE
SNE, I-SNE, UMAP
Well-known attraction/repulsion methods such as SNE [II] t-SNE 32 and
UMAP 19 define an affinity matrix with entries Pij in the ambient space, and another affinity matrix with entries dú in the embedded space. To define the embedding, a loss between the two affinity matrices minimized. Specifically, the loss function is DKL (pllq) == E.y Ptj Ptg Pty/93 in SNE and t-SNE, whereas UMAP is equivalent to adding DKL - pll - g) for unnormalized densities 22 While these methods preserve affinities, they do not preserve any types of distances in the embedding.

---

4 Heat-Geodesic Embedding
In this section, present our Heat Geodesic Embedding which is summarized in Alg. We start by introducing the heat-geodesic dissimilarity, then present a robust transformation, and a heuristic to choose the optimal diffusion time. Proofs not present in the main text given in AppendixA|
We consider the discrete case, where we have a set of points {2:) =: X in a high dimensional Euclidean space z; E Rd. From this point cloud, we want to define a map ø R + R that embeds the observation in lower dimensional space. An important property of our embedding is that we preserve manifold geodesic distances in low dimensional space.
Heat-geodesic dissimilarity Inspired by Varadhan's formula and the Harnack inequalities, we defined heat-geodesic dissimilarity based on heat diffusion graphs. From observations (data- points) in Rd, we define an undirected graph G, and compute its heat kemel H, = e e- tL, where L is the combinatorial symmetrically normalized graph Laplacian (the heat kernel is thus symmetric). Following the inequality (2) we can rearange the terms to isolate the geodesic distance, inspired by this observation, we define the following dissimilarity.
Definition 4.1. For a diffusion time > 0 and tunable parameter o > 0, we define the heat-geodesic dissimilarity between Fi, Ti X as
(x(,, 25)) == [-41 log(H))) o41 log Vi) v)
where H, is the heat kernel on the graph G, and (Vi)s := 2[(H))s; + (H)
Here the log is applied elementwise, and the term - 4t log HH) corresponds to the geodesic approximation when + 0 as in Varadhan's 's formula. In practice one uses a fixed diffusion time t > 0, so we add symmetric volume correction term as in the Harnack inequality, ensuring that de(z;, 2ty) is symmetric. From Sec. 2 we have he (2, x) = V(z, vt)-1, and we use the diagonal of H, approximate the inverse of the volume. With this volume correction term and 1, the dissimilarity is such that de(2,, x;) = 0 for all > O. When a = 0 or the manifold has uniform volume growth (as in the constant curvature setting) we show that the heat -eodeeeee dissimilarity is order preserving:
Proposition 4.2. When
0 or the manifold has uniform volume growth, i.e. (Ht)a: (Hz) as and the heat kernel is pointwise monotonically decreasing w.r.t. a norm |- in ambient space, we have for triples Zt, J,. 2x E X that [z, zz} > 12 zxl implies di (2;,,,,)) (e (k), i.e. the heat- geodesic dissimilarity is order preserving.
Proof. When
0 or the manifold has uniform volume growth we need only consider the -4t log He)k terms. The assumption of pointwise monotonicity of the heat kernel entails that [x; Zz > [r; zgl implies He( 2z) < HH( IK). We are able to conclude that 4t log (z,, 23) > 4t log H, (2,, mx) and thus dz(),,,))) > de(x;,
Denoising distances with triplet computations We note that both diffusion maps and PHATE compute a triplet distance between datapoints, i.c., rather than using the direct diffusion probability between datapoints, they use distance between corresponding rows of diffusion operator. In particular, difffusion maps uses Euclidean distance, and PHATE uses an M-divergence. Empirically, notice that this step acts as denoiser for distances. We formalize this observation in the following proposition. We note Dr the triplet distance. The triplet distance compares the distances relative to other points. Intuitively, this is a denoising step, since the effect of the noise is spread across the entire set of points. For a reference dissimilarity like the heat- geatdgee it is defined as Dr (r() := (I1,( de -12. For linear perturbations of the form de de mry) e, where e E R, the effect of e on linear combination between the neat-geodesic dissimilarity and its triplet distance where p E [0.
) is less severe than on (x 22). Our embedding is based on a
p)de + pDr,
Proposition 4.3. Denote the perturbed triplet distance by Dr 2z) = ||do(ss, -) de All
where de( x;) == de( x5)) + and (1 (I))) == (1( for k + j. Then the triplet distance Dr is robust to perturbations, i.e.,fr all e 0,
2
de [d((),,z)+ +
<
Dr Dr(2,,
d(()))2))
5

---

Optimal diffusion time Varadhan's formula suggests small value of diffusion time to approxi- mate geodesic distance on manifold. However, in the discrete data setting, geodesics are based on graph constructions, which in turn rely nearest neighbors. Thus, small can lead to disconnected graphs. Additionally, increasing can serve as a way denoising the kernel (which is often computed from noisy data) as it implements advantage adding noise tolerance. By computing a sequence of heat kernels (He)) and evaluating their entropy H (H)) ;= - ( log( Hi)ij. select & with the knee-point method [26] on the function H (H.). We show in Sec. 6.1 that heuristic for determining the diffusion time automatically leads to better overall results.
tow-pass filter over the eigenvalues, providing the additional
Weighted MDS The loss in MDS (eq is usually defined with uniform weights. Here, we optionally weight the loss by the heat kernel. In Sec. we will show how this modification relates our method to the embedding defined by SNE[D]] For Zi E X, we minimize (Ht) bi (de (1a, [[((;;) ((x)))l2)² This promotes geodesic preservation of local neighbors, since more weights are given to points with higher affinities.
Heat-geodesic embedding To define a lower dimensional embedding of point cloud X, we construct matrix from the heat-geodesic dissimilarity, and then MDS create the embedding. Our embedding defines a map @ that minimizes (d, (zt xyy) |l((x,) - (xx,)) lll) for all 24,, E X. Hence, it preserves the heat-geodesic dissimilarity as the loss decreases to zero. In Alg, we present the main steps of our algorithm using the heat-geodesic dissimilarity, detailed version is presented in Appendix A
Algorithm Heat Geodesic Embedding
1: Input: N d dataset matrix X, denoising parameter [0, 1]. Harnack regularization >> 0,
output dimension
2: Returns: N x k embedding matrix E.
3: He 4 Pk (L,
b Heat approximation Knee detection e.g. 20 log applied elementwise Triplet interpolation step
4: f Kneedle (H (.))).
5: D [-4 log( (H.);j - a4t log (V)),)
6: D (1 D) + pDT
7: Return E < MetricMDS(DD, -112,
Relation to other manifold learning methods
In this section, we elucidate theoretical connections between the Heat Geodesic Embedding and other manifold learning methods. We relate embeddings via the eigenvalues of H, or P with Laplacian eigenmaps and diffusion maps. We then present the relation between our methods and PHATE and SNE. We provide further analysis in the Appendix A In particular, introduce a new definition of kernel preserving embeddings; either via kernel-bace based distances (diffusion maps, PHATE) or via similarities (e.g. t-SNE, UMAP).
Diffusion maps with the heat kernel Diffusion maps define an embedding with the first eigenvectors (de); of P. while Laplacian eigenmaps [1 uses the eigenvectors (vi); of L. In the following, we recall the links between the two methods, and show that rescaled Laplacian eigenmaps preserves the diffusion distance with the heat kernel H,,
Lemma 5.1. Rescaling the Laplacian eigenmaps embedding with (e 12-211gh preserves the diffusion distance DM.
21 (b.
Relation to PHATE The potential distance in PHATE (Sec. 3 is defined by comparing the transition probabilities of two I-steps random walks initialized from different vertices. The transition matrix p mimics the heat propagation on graph. The heat- heat-geodesic dissimilarity provides a new interpretation of PHATE. In the following proposition, we show how the heat- nat--ge relates to the PHATE potential distance with a linear combination of t-steps random walks.

---

Proposition 5.2. The PHATE potential distance with the heat kernel РHH, can be expressed in terms of the heat- geodesic dissimilarity with o
РHн, (1/4t)² (dd(a..))) de de
and it is equivalent to multiscale random walk distance with kernel Σ mi (k)) where mk(k)= == the t /kkl.
Proof. We present a simplified version of the proof, more details are available in Appendix A For = O, we have dz(z;,z)) -4t -4llog' log (H) the relation between the PHATE potential and the heat-geodesic follows from the definition
PHH. (I),
Σ (- He( + log H((z,, (L/4t))lde(zz,.)) de(zs,
Using the heat kernel H₂ with the random walk Laplacian Ly to multiscale random walk kernel. We can write Lyw = SAS-1, where S = Q-1/2y Since P Ln Lyw we have p SI, S Interestingly, we can relate the eigenvalues of H; and P with the Poisson distribution. The probability mass function of a Poisson distribution with mean / is given by m,, relationship, we can express H, linear combination of ps weighted by the Poisson distribution. Indeed, substituting = f in yields
W corresponds
(k) == e /k:!. For 2 O, we have e ( =
Ek>o mz(k)) With this
00
H. = Se S S S
> (k)) - ( Σ my (k)
Remark 5.3. In the previous proposition, the same argument holds for the symmetric Laplacian and the affinity matrix A : Q21/ used in other methods such as diffusion maps (5]. This is valid since we can write Laym = Q--PP. 1/-/AAATTOO Q 1/2, and A = I, Lyg.-
Remark 5.4. This proposition shows that, as the denoising parameter P 1, Heat Geodesic Embedding interpolates to the PHATE embeddings with a weighted kernel Ex= my (k)
The heat- gea--geodesic method also relates SNE [1] and its variation using the Student distribution t-SNE [18]. In SNE, the similarity between points is encoded via transition probabilities Pij. The objective is learn an affinity measure g, that depends on the embedding distances |/y; 1/ [l2, such that it minimizes (pllg) Intuitively, points that have strong affinity in the ambient space, should also have a strong affinity the embedded space. Even though the heat-geodesic minimization is directly on the embedding distances, can show equivalent with SNE. Appendix [A] we provide additional comparisons between SNE and our method.
Relation SNE
Proposition 5.5. The Heat- Grat-Geedeescc embedding with a and squared distances minimization weighted by the heat kernel is equivalent to SNE with the heat kernel affinity in the ambient space, and a Gaussian kernel in the embedded space Qij =
- 1/ /
Results
In this section, show the versatility of our method, showcasing its performance in terms of clustering and preserving the structure of continuous manifolds. We compare the performance of Heat Geodesic Embedding with multiple state-of-the-arr baselines on synthetic datasets and real- world datasets. For all models, we perform sample splitting with 50/50 validation-test split. The validation and test each consists of repetitions with different random initializations. The hyper-pparae parameters are selected according to the performance on the validation set. We always report the results on the test set, along with the standard deviations computed over the five repetitions. We use the following methods in our experiments: our Heat Geodesic Embedding, diffusion maps [5], PHATE [22], shortest-path (used in Isomap 31) which estimates the geodesic distance by computing the shortest path between two nodes in graph built on the point clouds, 1-SNE (32). UMAP [19, and metric MDS with Fuclidean distance. Details about each of these methods, and results for different parameters (graph type, heat approximation, etc.) are given in Appendix C
7

---

Table 1: Pearson and Spearman correlation between the inferred and ground truth distance matrices on the Swiss roll and Tree datasets (higher is better). Best models on average are bolded.
Swiss roll
Tree
Distance
Pcarson
Spcarman
Pearson
Spcarman
Diffusion distance
0.476 + 0.226
0,478 = 0.138
0.656
0.054
0.653 + 0.057
PHATE potential
0.457 = 0.01
0.404 = 0.024
0.766 = 0.023
0.743 # 0.028
Shortest path
0.497 + 0.144
0.558 0.134
0.780 0.009
0.757 + 0.019
Fuclidean
0.365 0.006
0.413 0.005
0.735 0.014
0.704 + 0.033
leel-geodesic (ours)
0.702 1 0.086
0.700
0.073
0.822 0.008
0.807 1 0.016
6.1 Distance matrix comparison
We start by evaluating the ability of the different distances or dissimilarities to recover the ground truth distance matrix of point cloud. For this task, use the Swiss roll and Tree datasets, for which the ground truth geodesic distance is known. The Swiss roll dataset consists of data points sampled on a smooth manifold (see Fig. [. The Tree dataset is created by connecting multiple high-dimensional Brownian motions in a tree-shape structure. In Fig. T we present embeddings of both datasets. Our method recovers the underlying geometry, while other methods create artificial clusters or have too much denoising, Because we aim at a faithful relative distance between data points, we compare the methods according the Pearson and Spearman correlations of the estimated distance matrices with respect to ground truth. Results are displayed in Tab. [ We observe that Heat Geodesic Embedding typically outperforms previous methods in terms of the correlation with the ground truth distance matrix, confirming the theoretical guarantees provided in Sec. 4 & 2 Additional results such as computation time and correlation for different noise levels are available Appendix C
Optimal diffusion time In Section 4 we described heuristic to automatically choose the diffusion time based on the entropy of H, In Fig. 2 we show that the knee-point of + H (H:). corresponds to a high correlation with the ground distance, while yielding low approximation eиror of the distance matrix (mea- sured by the Frobenius norm of the difference between D and the ground truth).
1.-
2:h
205
07-
26-
6.2 Preservation of the inherent data structure
Eroww
Astwerik Tere
crucial evaluation criteria of manifold learning methods the ability capture the inherent structure of the data. For instance, clusters in the data should be visible in the resulting low dimensional representation. Similarly, when the dataset consists of samples taken at different time points, one expects be able characterize this temporal evolution in the low dimensional embedding (22). We thus compare the different embedding methods according to their ability retain clusters and temporal evolution of the data.
50
100
ffuiiion time
Figure 2: tion between estimated and ground truth distance matrices in function of the diffusion time t.
Evolution of the correla-
Isomap
Maps
Figure 3: Embeddings of 2000 differentiating cells from embryoid body [22] over 28 days. UMAP and t-SNE do not capture the continuous manifold representing the cells' evolution.
Identifying clusters. We use the PBMC dataset, the Swiss roll, the Tree dataset, MNIST g and COIL-20 (3 dataset. The PBMC dataset consists of single-cell gene expressions from 3000 individual peripheral blood mononuclear cells. Cells are naturally clustered by their cell type. For
8

---

PHATE
HeatSeo
HeatGeo 0.5
HeatGeo p
Call byso
co14+ Monoeytee
cels
c34 colk
Deed Meystayy-yat
Figure 4: Embeddings on PBMC using the triplet distance with the heat-geodesic for different regularization parameter p.
the Tree dataset, on the manifold according to mixture of Gaussians and use the mixture component as the ground truth cluster labels. The MNIST and COIL--20 datasets are clustered by digits objects but may not respect the manifold hypothesis. For each method, we run k-means on the two-dimensional embedding and compare the resulting cluster assignments with ground truth. Tab. [0] reports the results in terms of homogeneity and adjusted mutual information (aMI). Heat Geodesic Embedding competitive with PHATE and outperforms I-SNE and UMAP on all metries except on the MNIST and COIL-20 datasets. Yet, we show in Appendix C that all methods tend to perform equally well when the noise level increases. In Fig. we present the PBMC embeddings of PHATE and HeatGeo, showing that HeatGeo interpolates to PHATE for 1.
the branches as clusters. For the Swiss roll dataset, sample data points
Table 2: Clustering quality metrics for different methods. We report the homogeneity and the adjusted mutual information (aMI). Best models average are bolded (higher is better).
zoll
PSMC
MNITT Masi44)
COIL Masie4.)
Metted
Earreporatty
2b01
Horregmentt
aMI
aM1
-M11
Hereegeant
0.100
0.720
0.117
IL145-2
(156-2-0,015
0L528
0L7100-2:0:54
0.005
Tumrop
0.600 11 100
0.7:3 I.885
0.675 1.1144
6:081 - 1110
11:002 0:007
0:71 6.007
1.312 1-11
0.71 1-0.901
0.000
0.381
0:001 1..00
0.011
2:0003
0
0:236
0.206 0:0
0.550-2-0:042
Diffiction Maps
0.013 (CC))
0.065
0.311 1.110
0.354
00 (80)
0.006 6.001
8.006 0.997
1.556: 1-0.00
0.013 - 1.000
0.112 1 0.11
HearGen ovens)
01820 0..018
0.740
00014
0.761
0.051
0.766 0..51
037 6.009
0.785 1.0.
0.8% 01:00
08:01:0044
(.506) (.177)
Temporal data representation. For this task, we aim at representing data points from population observed at consecutive points in time. We single cell gene expression datasets collected across different time points, including the Embryoid Body (EB), IPSC (22). and two from the 2022 NeurlPS multimodal single-cell integration challenge (Cite & Multi). To quantitatively evaluate the quality of the continuous embeddings, we first embed the entire dataset and obfuscate all samples from a particular time point (e.g., 2). We then estimate the distribution of the missing time point by using displacement interpolation (35] between the adjacent time points (e.g.. = 1 and t 3). We report the Earth Mover Distance (EMD) between the predicted distribution and true distribution. A low EMD suggests that the obfuscated embeddings are naturally located between the previous and later time points, and that the generated embedding captures the temporal evolution of the data adequately. Results are presented in Tab. 3 Heat Geodesic Embedding outperforms other methods on the EB, Multi, and IPSC datasets and is competitive with other approaches on Cite. We show graphical depiction of the different embeddings for the embryoid (EB) dataset in Fig.
Table 3: EMD between a linear interpolation of two consecutive time points t - 1, f. + 1, and the time points f. Best models on average are bolded (lower is better).
Method
Cite
EB
Multi
IPSC
Euclidean
0.978 0.069
1.012 0.039
1,212 0.199
1.085 0.234
Isomap
0.978 = 0.105
0.993 - 0.062
1.299 - 0.307
1.026 0.253
Non-metric MDS
0.81 0.012
0.85 385+0.000 0.014
0.806 0.015
1.013 0.067
UMAP
0.791 = 0.045
0.942 = 0.053
1.418 0.042
0.866 + 0.058
I-SNE
0,905 = 0.034
0.964 = 0.032
1.208 = 0.087
1.006 + 0.026
PHATE
1.032 + 0.037
1.088 + 0.012
1.254 + 0.042
0.955 + 0.033
Diffusion Maps
0.989 - 0,080
0.965 = 0.077
1.227 - 0.086
0,821
- 0.039
HeatGeo (ours)
0.890 0.046
0.733 0.036
0.958 0.044
0.365 0.056

---

7 Conclusion and Limitations
The ability to visualize complex high-dimensional data in an interpretable and rigorous way is crucial tool of scientific discovery. In this work, we took step in that direction by proposing general framework for understanding diffusion-baar based dimensionality reduction methods through the lens of Riemannian geometry. This allowed us to define novel embedding based on the heat geodesic dissimilarity-a more direct measure of manifold distance. Theoretically, we showed that our methods brings greater versatility than previous approaches and can help gaining insight into popular manifold learning methods such as diffusion maps, PHATE, and SNE. Experimentally, we demonstrated that it also results in better geodesic distance preservation and excels both at clustering and preserving the structure of a continuous manifold. This contrasts with previous methods that are typically only effective at single of these tasks.
Despite the strong theoretical and empirical properties, our work presents some limitations. For instance, our method is based on similarity measure, which is relaxation of distance metric. Additionally, the Harnack equation suggests that our parameters for the volume correction could be tuned depending on the underlying manifold. We envision that further analysis of this regularization a fruitful direction for future work.
Acknowledgments and Disclosure of Funding
This research was enabled in part by compute resources provided by Mila (mila.quebec). It par- tially funded and supported by ESP Mérite [G.H.]. CIFAR AT Chair [G.W.]. NSERC Discovery grant 03267 [G. W.], NIH grants (1 F30A1157270000, R01HD100035, R01GM1 R01 RO1GM135929) GGW.S.SK.. NSF Career grant 2047856 [S.K.]. the Chan-Zuckerberg Initiative grants CZF2019- 182702 and CZF2019-002444 [. the Sloan Fellowship FG-2021- [S.K.]. and the Novo Nordisk grant GRI 12933 [S.K.]. The content provided here is soiely the responsibility of the authors and does not necessarily represent the official views of the funding agencies. The funders had no role in study design, data collection & analysis, decision to publish, or preparation of the manuscript.
References
[1] Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation, 15(6):1373-139g) 2003.
[2] Jan Niklas Böhm, Philipp Berens, and Dmitry Kobak. Attraction-repulsion spectrum in neighbor embeddings. Journal of Machine Learning Research, 23(95):1-32, 2022.
[3] Erica I. Busch, Jessie Huang, Andrew Benz, Tom Wallenstein, Guillaume Lajoie, Guy Wolf, Smita Krishnaswamy, and Nicholas BrownBr Multi-view manifold learning human brain-state trajectories. Nature Computational Science, pages 1-14, 2023.
[4] J Douglas Carroll and Phipps Arabie. Multidimensional scaling. Measurement, judgment and decision making, pages 179-250, 1998.
[5] Ronald R. Coifman and Stéphane Lafon. Diffusion maps. Applied and Computational Harmonic Analysis, 21(1):5-30 2006.
[6] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. Geodesics in heat: A new approach to computing distance based on heat flow. ACM Transactions Graphics (TOG), 1111, 2013.
[7] Michaël Defferrard, Lionel Martin, Rodrigo Pena, and Nathanaël Perraudin. Pygsp: Graph signal processing in python.
[8] Li Deng. The mnist database of handwritten digit images for machine learning research. IEEE Signal Processing Magazine, 29(6): 141-142, 2012.
[9] Laleh Haghverdi, Maren Büttner, F. Alexander Wolf, Florian Buettner, and Fabian J. Theis. Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods, 13(10):845-848, 2016.
10

---

[10] Matthieu Heitz, Nicolas Bonneel, David Coeurjolly, Marco Cuturi, and Gabriel Peyré. Ground metric learning on graphs. Journal of Mathematical Imaging and Vision, 63:89-107, 2021.
[11] Geoffrey E Hinton and Sam Roweis. Stochastic neighbor embedding. Advances neural information processing systems, 15, 2002.
[12] Michael Hout, Megan H Papesh, and Stephen D Goldinger. Multidimensional scaling, Wiley Interdisciplinary Reviews: Cognitive Science, 4(1):93-103, 2013.
[13] Shih-Gu Huang, Iwwo Lyu, Anqi Qiu, and Moo K Chung. Fast polynomial approximation of heat kernel convolution on manifolds and its application to brain sulcal and gyral graph pattern analysis. IEEE transactions on medical imaging, 39(6):2211-22) 2020.
[14] Guillaume Huguet, Alexander Tong, Bastian Rieck, Jessie Huang, Manik Kuchroo, Matthew Hir,, Guy Wolf, and Smita Krishnaswamy. iime-inhomogeneou diffusion geometry and topology. arXiv preprint arXiv:2203. 14860, 2022.
[15] Guillaume Huguet, Alexander Tong, María Ramos Zapatero, Guy Wolf, and Smita Krish- naswamy. Geodesic Sinkhorn: optimal transport for high- digh-dimensioal datasets. arXiv preprint arXiv:2211.00800 2022.
[16] Joseph Kruskal. Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. Psychometrika, 29(1):1-27, 1964.
[17] Peter Li and Shing Tung Yau. On the parabolic kernel of the Schrödinger operator. Acta Mathematica, 156(none):11 1533 201, 1986.
[18] George C. Linderman and Stefan Steinerberger. Clustering with t-SNE, provably. arXiv:1706.0223 [cs, stat], 2017.
[19] Leland Mclnnes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:/802.03426,, 2018.
[20] Stanislav A Molchanov. Diffusion processes and riemannian geometry. Russian Mathematical Surveys, 30(1):1, 1975.
[21] Kevin R. Moon, Jay S. Stanley, Daniel Burkhardt, David van Dijk, Guy Wolf, and Smita Krishnaswamy, Manifold learning- based methods for analyzing single-cell RNA-sequencing data. Current Opinion in Systems Biology, 7:36-46, 2018.
[22] Kevin R. Moon, David van Dijk, Zheng Wang, Scott Gigante, Daniel B. Burkhardt, William S. Chen, Kristina Yim, Antonia van den Elzen, Matthew J. Hir, Ronald R. Coifman, Natalia B. Ivanova, Guy Wolf, and Smita Krishnaswamy. Visualizing structure and transitions in high- dimensional biological data. Nat Biotechnol, 37(12): 182-11422 2019.
[23] Samer A. Nene, Shree Nayar, and Hiroshi Murase. Columbia object image library (coil- 20). Technical Report February 1996.
Department of Computer Science, Columbia University,
[24] Adam Nowak, Peter Sjögren, and Tomasz 7. Szarek. Sharp estimates of the spherical heat kernel. Journal de Mathématiques Pures et Appliquées, 129:23-33, 2019.
[25] Laurent Saloff-Coste The heat kernel and its estimates. Probabilistic approach to geometry, 57:405-436, 2010.
[26] Ville Satopaa, Jeannie Albrecht, David Irwin, and Barath Raghavan. Finding a" kneedle" in a haystack: Detecting knee points in system behavior. In 2011 31st international conference distributed computing systems workshops, pages 166-171. IEEE, 2011.
[27] Nicholas Sharp, Yousuf Soliman, and Keenan Crane. The vector heat method. ACM Tronsactions on Graphics (TOG), 38(3):1-19, 2019.
[28] Justin Solomon, Fernando De Goes, Gabriel Peyré, Marco Cuturi, Adrian Butscher, Andy Nguyen, Tao Du, and Leonidas Guibas. Convolutional wasserstein distances: Efficient optimal transportation geometric domains. ACM Transactions on Graphics (ToG), 34(4):1-11, 2015.
11

---

[29] Jian Sun, Maks Ovsjanikov, and Leonidas Guibas. A concise and provably informative multi- scale signature based on heat diffusion. In Computer graphics forum, volume 28, pages 1383-1392. Wiley Online Library, 2009.
[30] Yoshio Takane, Forrest W Young, and Jan De Leeuw. Nonmetric individual differences multidimensional scaling: An alternating least squares method with optimal scaling features. Psychometrika, 42:7-67, 1977.
[31] Joshua B Tenenbaum, Vin de Silva, and John C Langford. A global geometric framework for nonlinear dimensionality reduction. science, 290(5500):231-22222 2000.
[32] Laurens Van der Maaten and Geoffrey Hinton. Visualizing data using t-sne. Journal of machine learning research, 9(11), 2008.
[33] David Van Dijk, Roshan Sharma, Juozas Nainys, Kristina Yim, Pooja Kathail, Ambrose Carr, Cassandra Burdziak, Kevin R Moon, Christine I. Chaffer, Diwakar Pattabiraman, et al. Recovering gene interactions from single-cell data using data diffusion. Cell, 174(3)716-729, 2018.
[34] Sathamangalam Srinivasa Varadhan. On the behavior of the fundamental solution of the heat equation with variable coefficients. Communications on Pure and Applied Mathematics, 20(2): 431-455, 1967.
[35] Cédric Villani. Displacement interpolation. Optimal Transport: Old and New, pages 113-162, 2009.
[36] F Alexander Wolf, Philipp Angerer, and Fabian J Theis. Scanpy: large-scale single-cell gene expression data analysis. Genome biology, 19:1-5, 2018
[37] Jonathan Zheng, Samraat Pawar, and Dan FM Goodman. Graph drawing by stochastic gradient descent. IEEE transactions visualization and computer graphics, 250)):738-2748, 2018.
12

---

Appendix
A Theory and algorithm details
A.1 Kernel preserving embeddings
In this section, we attempt to create a generalized framework for dimensionality reduction methods. These methods often have been viewed as disparate or competing but here we show that many of them are related to one another given the right template for methodology comparison. In order to do this, we introduce a general definition suited for With this definition, we can cast many dimensionality reduction methods within the same framework, and easily compare them. We recall that the observations in the ambient space are denoted æ. and those in the embedded space are denoted y. The definition relies on kernel functions H;"" H2 defined respectively on the ambient and embedded spaces and on transformations T TM applied to the kernels. We recall that divergence f IR x IK - R+ is such that [(a b) 0 if and only if a = b and (a, a + 8) a positive semi-definite quadratic form for infinitesimal
dimensionality reduction methods.
Definition A.1. We define kernel features preserving embedding as an embedding which minimizes loss L between a transformation 2 of the ambient space kernel Hü and its embedded space counterpart
L : f(T Th (H)),,
(6)
where f is any C2 divergence on Rzo-
Example 1. We formulate MDS kernel feature-preserving embedding, Suppose we want to preserve the Euclidean distance, we have (Z = ||< xzl12. H² (y)) = l - y/ 35ll2.
fa, = ||a bll2, and T T
In the following, we present popular dimensionality reduction methods that are kernel features preserving embeddings. With this definition, we can distinguish between methods that a preserve kernel via affinities or distances. For the methods considered in this work, H. is an affinity kernel, but its construction varies from one method to another. In PHATE and Diffusion maps, H₂" is random walk P, while in Heat Geodesic Embedding we use the heat kernel H. t-SNE defines H. as symmetrized random walk matrix from Gaussian kernel, while UMAP uses an unnormalized version. Methods such as PHATE and diffusion maps define new distance matrix from kernel in the ambient space and preserve these distances in the embedded space. Other methods like t-SNE and UMAP define similarities from a kernel and aim to preserve these similarities in the ambient space and embedded space via antropy-basce loss. We note the Kullback-Leibler divergence DKL (a, b) = Eij aiş log(a.;/biz].
Proposition A.2. The embeddings methods HeatGeo, PHATE, Diffusion Maps, SNE, I-SNE, and UMAP are kernel cerature-presenvee embeddings.
Proof. We assume that the affinity kernel in the ambient space H! is given, to complete the proof we need to define f, H,,," T=,T TT for all methods.
We start with the distance preserving embeddings; HeatGeo, PHATE, and Diffusion Maps. For these methods, the kernel in the embed space is simply H (y4. U3y) = 11y, - y (2. without transformation, i.e T T. Since they preserve distance, the loss is [(T" (H;) T (H²)) = !| H ||2-
In the Heat Geodesic Embedding we apply a transformation on H H, to define dissimilarity, hence T (HT) = -tlog H. (for a = 0), where log is applied elementwise.
In PHATE, the potential distance is equivalent (I (HFF) = || -| log( + log(H)())) In Diffusion Maps, the diffusion distance is (IT(H)) (HT)) = [1(11)); (HrF),l2.
SNE, t-SNE, and UMAP preserve affinities from a kernel. For these three methods, the loss is a divergence between distributions, namely f = DKL- They vary by defining different affinity kernel and transformation in the embedded space. SNE uses the unnormalized kernel H.¹ (vt, v,)
exp(=(1/t)llu 2)) with T = T I. Whereas, t-SNE uses (H1)) (1 + |/ 1
13

---

Table 4: Overview of kernel preserving methods.
Method
H (y:. m₁)
T (H;;")
T (H
f
PHATE
M/ Vy 2
log H log( (Hg(H)))
H
Heat Geodesic
yi Uy 2
-log H3 kij
Hy
Diffusion Maps
Iy 3/ 2
[1HH;):
(H
H
1-112
SNE
exp( - ( I my 1))
H
H.
DKL
I-SNE
(1 + 1& 1/
H
DKL
UMAP
(1 + 1& y/ )
H
(8P)
DKL
(1-(H))6)
and 7 = 7 = 1. UMAP define a pointwise transformation in the embedded space with (1T))
(1 + 1/y/ - y, ( (H,))) = (1 - (HT) and T I.
We summarize the choice of kernels and functions in Tab. 4
A.2 Proofs
Proposition 4.3. Denote the perturbed triplet distance by Dr Fz) (1d((s,,)) de(syy l2
where de(x;; k5) == de(x;; x5) + e and de (2,, (x)) == (e, 2k.) for k # j. Then the triplet distance DT is robust to perturbations, ie., for all > 0.
DT(
de (1) )
Dr(Z(,Iz)
Proof of Proposition 4.3 The effect of the noise on the square distance is (de(x;; Zz)
e)² /d(x;, 2,)² 1 (2ed, 24) + 2)),,,, Denoting the perturbed triplet distance by Dr. we have
4ed(s;, mz) + 202 Extis (de( - xx)x 2 + 2(d + =1+ Dr(E,,I,) Dr(I;,2,)²
o,) Dr(r,,2,)²
and we have
4ed( 2;) x 2
2ed (x;,, +
DT(11,
de(
For >> 0, this gives
4d (31, x.y) 2dd 2de 2r)
(x:;
Di(E,, x,,² - 2dl (1;,,
For < have
4d Adi 2dh (d ( D. De xz2)² ( (I,) x,)²
(I.
Thus E E (-00, (x, (0, 00). As we require the perturbation factor e << de
hence choose (0,)0).
Lemma 5.1. Rescaling the Laplacian eigenmaps embedding with v; + (2-2 2t tbk.)) preserves the diffusion distance DMH.
Proof of Lemma 5.7 Since the eigendecompatision of H, form orthonormal basis of T and
since its first eigenvector is constant, we can write the diffusion distance ||,, 5, H, Σx20 E 22 (U's -xx,)² Ex21 -11 (b' In particular, this defines the dimensional embedding + (t (x)...., tAx (b)))(()))
14

---

Proposition 5.2. The PHATE potential distance with the heat kernel РHH, can be expressed in terms of the heat- geodesic dissimilarity with
РHH, (1/4t))[dd - de (2z. 1l1
and is equivalent to multiscale random walk distance with kernel Ekso mi (i) where mi == t=tet
Proof of Proposition For O, we have de (Zz, zy) 4t log (He)s. the relation between the PHATE potential and the heat-geodesic follows from the definition
РHН. = Σ - log He + log Hi Hi(zz, xx))
k
(L/44)
Using the heat kerne H, with the random walk Laplacian Lyw a multiscale random walk kernel. Recall that
corresponds write Lew terms of the symmetric Laplacian L = Q-1/2 L QL/2 meaning that the two matrices are similar, hence admit the same eigenvalues A. We also know that Lg is diagonalizable, since
write L, Q-11-QQQ----------
Q 0-1/24A Q 11/2 In particular, we have Lyw SAS where S == 11 The random walk matrix can be written as P I, Rrw, hence eigenvalues are (In A), and we can write p === S - A Similarly, the heat kernel with the random walk Laplacian can be written as H. = Se-tAS s-1 Interestingly, we can relate the eigenvalues of H, and P with the Poisson distribution. Note the probability mass function of a Poisson as m,(k) == e //l, for 2 0, we have
(tu)*
((1--) =
Σ me(k)y²
k!
$20
k20
We note that (7) is the probability generating function of a Poisson distribution with parameter f, i.e. J, where X Poisson(t). With this relationship, we can express H, as linear combination of P weighted by the Poisson distribution. Indeed, substituting À = in 77 links the eigenvalues of H, and P. We write heat kernel as a linear combination of random walks weighted by the Poisson distribution, we have
H. Se Se = S E mz (L) - -A) Σ mi (k) k=
k=0
Proposition S.5. The Heat-Geaue Geea-GGedeesic embedding with 0 and squared distances minimization weighted by the heat kernel equivalent to SNE with the heat kernel affinity the ambient space, and a Gaussian kernel in the embedded space qı exp(-lly.
/4f).
Proof of Proposition 5.5 The squared MDS weighted by the heat kernel corresponds to
Σ
by (Ti, 23 (d) - 1ly. yÿ 1) = Σ
by (x my) ( (- log he (z) - llys yز 1)
1
15
Σ
he (xs, v, )t² (log he (x), log exp(-lla) 11112//27 /1/2
=
there exists an embedding that attain zero loss, then is the same as Eng he (1),, hel,, 2r) engexp(-llee 13 112/1)) DRr. (hx [lg).
A.3 Algorithm details
We present a detailed version of the Heat Geodesic Embedding algorithm in Alg
For the knee-point detection use the Kneedle algorithm [26]. It identifies a knee-point point where the curvature decreases maximally between points (using finite differences). We summarize the four main steps of the algorithm for function fæz), and we refer to (26) for additional details.
15

---

Algorithm 2 Heat Geodesic Embedding
Input: N d dataset matrix X, denoising parameter p [0, 1J, Harnack regularization 0,
output dimension k.
2: Returns: x embedding matrix E.
3: & I. Calculate Heat Operator H
4: if is "auto" then
5: i + Kneedle )
Knee detection e.g.
26
6: W kernel(X)
7: L <- Q 3-1 W
8: if Exact then
9: H, V
10: else
11:
H PK (L,t)
12: b Calculate Pairwise Distances D
13:
D --44llog log H,
log is applied elementwise
14:
D (1-1- p) pDT
Triplet interpolation step
Return E < MetricMDS(D, || 112,,
1. Smoothing with a spline to preserve the shape of the function.
2. Normalize the values, so algorithm does not depend on the magnitude of the observations.
3. Computing the set of finite differences for x and y (x), e.g. y4, : f(x;) I₁.
4. Evaluating local maxima of the difference curve yd., and select the knee-point using threshold based on the average difference between consecutive æ.
B Experiments and datasets details
Our experiments compare approach with multiple state-of-the-art baselines for synthetic datasets (for which the true geodesic distance known) and real-world datasets. For all models, we perform sample splitting with 50/50 validation-test split. The validation and test sets each consist of repetitions with different random initializations. The hyper-parameters are selected according to the performance the validation set. We always report the results the test along with the standard deviations computed over the five repetitions. We use the following state-of-the-an methods in our experiments: Heat Geodesic Embedding, diffusion maps [5), PHATE [22]. Heat-PHATE (a variation of PHATE using the Heat Kernel), Rand-Geo (a variation of Heat Geodesic Embedding where the random walk kernel), Shortest-path which estimates the geodesic distance by computing the shortest path between two nodes graph built on the point clouds, 1-SNE and UM
B.1 Datasets
We consider two synthetic datasets, the well known Swiss roll and the tree datasets. The exact geodesic distance can be computed for these datasets. We additionally consider real-world datasets: PBMC, IPSC [22), [22), and two from the from the 2022 NeurIPS multimodal single-cell integration challenge
B.1.1
Swiss Roll
The Swiss roll dataset consists data points samples on a smooth manifold inspired by shape famous alpine pastry. In its simplest form, it is 2-dimensional surface embedded in R given by
I L. cos
h
= t. sin(t)
16

---

where t [To.Ti] and h E [0, W]. In our experiments we used Tú = = T1 nT, and W = 5. We use two sampling mechanisms for generating the data points ; uniformly and clustered. In the first, we sample points uniformly at random in the (To. Ti x [0, W] plane. In the second, we sample according to a mixture of isotropic multivariate Gaussian distributions in the same plane with equal weights, means [(7, W/2), (12, W/2)], and standard deviations [1, 1]. In the clustered case, data samples are given a label according to the Gaussian mixture component from which they were
sampled.
We consider variations the Swiss roll by projecting the data samples in higher dimension using random rotation matrix sampled from the Haar distribution. We use three different ambient dimensions: 3, 10, and 50.
Finally, we add isotropic Gaussian noise to the data points in the ambient space with standard deviation
B.1.2 Tree
The tree dataset is created by generating K' branches from a D-dimensional Brownian motion that eventually glued together. Each branch is sampled from a multidimensional Brownian motion dX = 2dW(t (t) at times t = 0, 1, 2, .... L - for k E [K]. The first branch is taken as the main branch and the remaining branches are glued to the main branch by setting Xx Xx Xoligl where 4 is random index of the main branch vector. The total number of samples is thus L K
In experiments, used 500, K and D 5, 10 (i.e. two versions with different dimensions of the ambient space).
B.2 Evaluation Metrics
We compare the performance of the different methods according several metrics. For synthetic datasets, where ground truth geodesic distance is available, we directly compare the estimated distance matrices and ground truth geodesic distance matrices. For real-world datasets, clustering quality and continuous interpolation as evaluation metrics.
B.2.1 Distance matrix evaluation
The following methods use an explicit distance matrix: diffusion maps, Heat Geodesic Embedding. Heat-Phate, Phate, Rand-Geo and Shortest Path. For these methods, we compare their ability their ability recover the ground truth distance matrix several metrics. Letting D and Ó the ground truth and inferred distance matrices respectively, and N the number points in the dataset, we use following metrics.
Pearson p We compute the average Pearson correlation between the rows of the distance matrices, N Em E where Tx.y is the Pearson correlation coefficient between vectors and D; stands for the i-th row of D.
Spearman P We compute the average Spearman correlation between the rows of the distance matrices, Ei D.. where T2, is the Spearman correlation coefficient between vectors æ and y. D, stands for the i-th row of D.
Frobenius Norm We use ||D - DlLF, where ||A||= = v E EA
Maximum Norm We use || DIloo where ||A] = max; |Ai5l
B.2.2
Embedding evaluation
Some methods produce tow-dimensional embeddings without using an explicit distance matrix for the data points. This is the case for UMAP and t-SNE. To compare against these methods, we use the distance matrix obtained by considering Fuclidean distance between the dow-dimensional embeddings. We used 2-dimensional embeddings in our experiments. For diffusion maps, obtain
17

---

these embeddings by using the first two eigenvectors of the diffusion operator only. For Heat Geodesic Embeddingg Heat-PHATE, PHATE, Rand-GEO and Shortest Path, we use multddmensional scaling (MDS) on the originally inferred distance matrix.
Clustering We evaluate the ability of Heat Geodesic Embedding create meaningful embeddings when clusters are present in the data. To this end, we run a k-means clustering on the two dimensional embeddings obtained with each method and compare them against the ground truth labels. For the Tree dataset, we use the branches as clusters. For the Swiss roll dataset, we sample data points on the manifold according to a mixture of Gaussians and use the mixture component as the ground truth cluster label.
Interpolation To quantitatively evaluate the quality of the continuous embeddings, first embed the entire dataset and obfuscate all samples from a particular time point (e.g.. = 2). We then
estimate the distribution of the missing time point using displacement interpolation B5 between
the adjacent time points (e.g.. t and = 3). We report the Earth Mover Distance (EMD) between the predicted distribution and true distribution. A low EMD suggests that the obfuscated embeddings are naturally located between the previous and later time points, and that the generated embedding captures the temporal evolution of the data adequately.
B.3 Hyperparameters
In Table 5
we report the values of hyperparameters used compute the different embeddings.
Hyperparameter
Description
Values
Heat Geodesic Embedding
k
Number of neighbours k-NN graph
5,10,15
order
order of the approximation
30
t
Diffusion time
0.1,1,10,50,aut
Approximation method
Approximation method for Heat Kernel
Buler, Chebyshev Combinatorial 0.0.2550.550075.11155
Laplacian
Type of laplacian
Harnack p
Haraack Regularization
PHATE
n-PCA
Number of PCA components
50,100
Diffusion time
1,5,10,20,auto
k
Number of neighbours
10
Diffusion Maps
k
Number of neighbours in k-NN graph Diffusion time
5,10,15
t
1,5,10,20
Shortest Path
k
Number of neighbours in k-NN graph
5,10,15
UMAP
k min-dist
Number of neighbours
5,10,15
Minimum distance
0.1,0.5.0.99
t-SNE
p early exageration
Perplexity
10,30,100
Early exageration parameter
12
Table 5: Hyperparameters used in our experiments
B.4 Hardware
The experiments were performed on a compute node with 16 Intel Xeon Platinum 8358 Processors and 64GB RAM.
18

---

C Additional results
C.1 HeatGeo weighted
Following Sec. [ we know that weighting the MDS loss by the heat kernel corresponds to a specific show the
parametrization of SNE, and thus promote the identification of cluster. In Fig- embeddings of four Gaussian distributions in 10 dimensions (top), and the PBMC dataset (bottom). The reference embedding using t-SNE, it models it also minimizes the KI. between the ambient and embedded distributions. We see that HeatCieo weighted form cluster that are shaped like Gaussian. This expected Prop. 5.5 indicates that this equivalent to minimizing the DKL between the heat kernel and a Gaussian affinity kernel.
HealGeo
HeatGeo Weighted
10D
Gaussians
PBMC
Figure 5: Embeddings of four Gaussian distributions in 10 dimensions (top), and the PBMC dataset (bottom). HeatGeo with weight is equivalent to minimizing the DxL. between the heat kernel and Gaussian affinity kernel, hence produces clusters shaped similar to Gaussian.
C.2 Truncated distance
In Fig 6, we discretize the interval [0, 51] in 51 nodes, and we compute the heat-geodesic distance of the midpoint with respect the other points, effectively approximating the Euclidean distance. Using Chebyshev polynomials of degree of 20, we see that the impact of the truncation is greater as the diffusion time increases. The backward Euler methods does not result in truncated distance.
Squared Heat-Geodesic distance 22
150
1500
Exact
Chebyshev
100
1000
Backward
4000
Harnack Reg
50
500
25
75
1.00
23
75
120
25
50
75
100
Figure 6: Approximation of the squared Euclidean distance with the Hcat-geodesic for the exact computation, Backward Buler approximation, and Chebyshev polynomials. For larger diffusion time, the Chebyshev approximation results in a thresholded distance, The Harnack regularization unsures (1 x) = O.
19

---

order
ooder
order
odde 30
12
order
order
onter
indo
undo
CIERY
crter 2>
onter 21
order 24
order 35
Figure 7: Impact of the Checbyshev approximation order on the embedding of HeatGeo for the PBMC dataset.
C.3 Harnack inequality
For complete Riemannian manifolds that satisfy the parabolic Harnack inequality (PHI) we have he(2,8) v) = (a, Vt) dd((,))) so that -1 log h (x v) = I log V(x, vt) + dd(r, y) 25
h1(2, x)) = V(x, VF) V(x, vT) -)((x,2)-1
(8)
(9)
We then have,
(xx v) ~ t ho v) t log Var, vt) 8 (z. 4) ~ -t log he(2, y) t llog he xx)
(x, y) ~ log (z. 8) + t tlog
C.3.1 Case studies for specific manifolds
The circle- SS S We now show that our expression for the Heat Geodesic Embbedding-distance is monotonically increasing with respect to the ground truth geodesic distance E IR+ for a fixed diffusion time t and for any Harnack regularization in Sj.
20

---

Our expression for the Heat Geodesic fimeedding-distamme is
V -4t log(h( 4/ 4/ log(h+(0))
As the square-root is monotonic, and 4t log h4 (0) is constant with respect to d, we need to show that f(d) (d) is monotonically increasing. For 1, we have
he(d) Σ MEZ
(2+22m)
VrT
(d+2cm))
As log is monotonically increasing, it suffices to show that Σ me decreasing, which is the case as for any d' > d, Vm E Z, we have
is monotonically
(d) 2xmm)
2
48
In general, one can see that (1) the heat kernel depending only on the geodesic distance and (2) the heat kernel being monotonically decreasing with respect to the geodesic distance sufficient conditions for preserving ordering of pair-wise distances with Heat Geodesic Embedding.
The sphere- Sn The above result can be applied the that the heat kernel on manifold of constant curvatures a function of the the geodesic distance (d) and time only. For Sn the heat kernel is given by
sphere S,, It is known
8
21 cc=-(cco(dd))
h, v)
1((+))-2
(cos(d))
n n- 2
1=0
with / the regularized incomplete beta function and C the Gegenbauer polynomials.
Furthermore, Nowak et al. (24] showed that the heat kernel of the sphere is monotonically decreasing.
The distance inferred from Heat Geodesic Embedding thus preserves ordering of the pair-wise distances.
Euclidean (IR3) For the euclidean space, we have for the volume of Vi-eeodssic ball and for the heat kernel:
4 T///
V V
he(x, (4
Recalling Harnack inequality.
dix.p? ug
d(m.v)?
S he(2,y) 5
V(x, v
V(x, vt)
With C2 Cy 4, we have
21

---

61 V(z, vt)
03
S (4Tt)3/2 $ ((z, Vt)
In this case, the bound can be made tight, by setting
Var, vl)
C1 C3
(4Tt)3/2
TT3/2
(Art)3//2
3v4T 6V
we recover the exact geodesic distance.
C.4
Quantitative results
C.4.1
Distance matrix evaluation
We report the performance of the different methods in terms of the ground truth geodesic matrix reconstruction in Table. for the Swiss roll dataset and Table. 7 for the Tree dataset.
C.4.2 Distance matrix evaluation via two-dimensional embeddings
We report the performance of the different methods in terms of the ground truth geodesic matrix reconstruction in Table for the Swiss roll dataset and in Table for the Tree dataset.
C.4.3 Clustering quality evaluation
On Tables 10, we report the performance on clustering quality for the synthetic datasets with different noise level.
C.5
Impact of the different hyperparameters
We investigate the impact of the different hyperparameters on the quality of the embeddings. In Figure R we show the embeddings of HeatGeo for different values of diffusion time, number of neighbours, order, and Harnack regularization.
In Figures 9 10 and 12 we show the impact different hyperparameters on the Pearson correlation between the estimated distance matrix and ground truth distance matrix for different methods on the Swiss roll dataset.
C.6 Graph construction
We compare the embeddings of the heat- geat-geodesic distance for different graph construction. Through- out the paper we used the graph construction from PHATE [22). In the following we present additional results depending on the choice of kernel construct the graph. Specifically, neighbor (kNN) graph implemented in 77), the graph from UMAP [19], and the implementation in the package Scanpy (36) for single-cell analysis. In figure, we present the embeddings 2500 points of tree with five branches in 10 dimensions, where the observations are perturbed with a standard Gaussian noise. All methods used five nearest neighbors and a diffusion time 20. Figure [3 we show the evolution of the Pearson correlation between estimated and ground truth distance matrices for the 10-dimensional Swiss roll dataset for various graph constructions, We note that the results are stable across different graph construction strategies.
use simple nearest
22

---

data
Noise level
Method
PearsonR
SpearmanR
Norm Fro N2
Norm inf N2
Swiss roll
0.1
Diffusion Map
0.974 + 0.01
0.983 0.007
0,018 +
0,026 + 0.0
Swiss roll
0.1
Hcat-Geo
0.992 + 0.003
0.995 0.002
0.002 + 0.0
0.003 = 0.0
Swiss roll
0.1
Heal-PHATE
0.99 0.002
0.997 = 0.001
0.079 0.002
0.1 0.003
Swiss roll
0.1
PHATE
0.621 + 0.006
0.58 0.01
0.022 + 0.0
0.026 + 0.0
Swiss roll
0.1
Rand-Geo
0.956 +: 0.003
0.993 0.001
0.009 + 0.0
0,012 + 0.0
Swiss roll
0.1
Shortest Path
1.0 + 0.0
1.0 0.0
0.0 0.0
0.001 + 0.0
Swiss roll
Euclidean
0.379 + 0.003
0.424 0.003
0.014 0.0
0.018 + U.O
Swiss roll
Diffusion Map
0.982 + 0.003
0.987 0.002
0.018 0.0
0.026 + 0.0
Swiss roll
0.5
Heat-Geo
0.994 + 0.002
0.996 0.001
0.002 0.0
0.004 + U.O
Swiss roll
0.5
Hcat-PHATE
0.993 + 0.001
0.998 = 0.0
0.064 + 0.001
0.083 = 0.002
Swiss roll
0.5
PHATE
0.649 + 0.007
0.615 0.006
0.023 + 0.0
0.028 + 0.0
Swiss roll
Rand-Goo
0.969 +: 0.002
0.995 0.001
0.009
0.011 = 0.0
Swiss roll
0.5
Shortest Path
0.999 + 0.0
0.999 0.0
0.001
0.002 + 0.0
Swiss roll
0.5
Euclidean
0.376 +: 0.004
0.422 0.004
0.0
0,018 = 0.0
Swiss roll
1.0
Diffusion Map
0.476 +: 0.226
0.478 = 0.138
0.018
0.026 + 0.0
Swiss roll
1.0
Heat-Geo
0.702 0.086
0.073
0.01
0.012 + 0.0 0.002
Swiss roll
Hcat-PHATE
0.623 +: 0.144
0.633 0.114
0.019 = 0.001
Swiss roll
1.0
PHATE
0.457 + 0.01
0.404 0.024
0.024
0,028 + 0.0
Swiss roll
1.0
Rand-Geo
0.521 +: 0.042
0.608 0.025
0.0
0.014 0.0
Swiss roll
1.0
Shurlest Path
0.497 : 0.144
0.558 0.134
0.011 0.001
0.015 + 0.002
Swiss roll
1.0
Euclidean
0.365 + 0.006
0.413 0.005
0.013 + 0.0
0.019 + 0.001
Swiss roll high
0.1
Diffusion Mup
0.98 0.003
0.986 0.001
0.018 0.0
0.026 + 0.0
Swiss roll high
0.1
Hcat-Geo
0.992 +: 0.003
0.996 0.002
0.002 + 0.0
0.003 + 0.0
Swiss roll high
0.1
Heat-PHATE
0.991 + 0.002
0.997 0.001
0.079 + 0.002
0.101 + 0.004
Swiss roll high
PHATE
0.625 + 0.013
0.582 0.017
1.022
0.026 = 0.0
Swiss roll high
0.1
Rand- Geo
0.956 + 0.002
0.993 0.001
0.009 + 0.0
0.012 + 0.0
Swiss roll high
0.1
Shortest Path
0.0
1.0
0.0
0.002 1 0.0
Swiss roll high
0.1
Euclidean
0.379 + 0.002
0.424 = 0.002
0.014 0.0
0.018 + 0.0
Swiss roll high
0.5
Diffusion Map
0.98 0.98 0.002
0.985 = 0.002
0.0
0.026 + 0.0
Swiss roll high
Hcat-Goo
0.997 + 0.001
0.997 0.0
0.005 0.0
0.007 0.0
Swiss roll high
0.5
Heat-PHATE
0.995 + 0.0
0.997 0.0
0.041 00111000
0.054 = 0.002
Swiss roll high
0.5
PHATE
0.717 0.004
0.707 = 0.005
0.0
0.034 = 0.001
Swiss roll high
0.5
Rand-Geo
0.984
0,996 0.0
0.008 0.0
0.01 + 0.0
Swiss roll high
Shortest Path
0.999
0.998 0.0
0.006 1
0.009 + 0.0
Swiss roll high
0.5
Euclidean
0.369 +: 0.003
0.421 0.003
0.013 + 0.0
0,018 + 0.0
Swiss roll high
1.0
Diffusion Map
0.555 +: 0.155
0.526 0.081
0.018
0,026 + 0.0
Swiss roll high
1.0
Heat-Geo
0.705 1: 0.065
0.695 0.052
0.011 + 0.0
0.012 : 0.0
Swiss roll high
1.0
Heal-PHATE
0.33::01066
0.625 0.074
0.011 1011:000 0.001
0.014 + 0.002
Swiss roll high
1.0
PHATE
0.473 + 0.026
0.419 0.024
0.027
0.0
0.039 + 0.001
Swiss roll high
1.0
Rand-Geo
0.563 +: 0.05
0.644 0.033
0.01 0.0
0.012 + 0.0
Swiss roll high
1.0
Shortest Path
0.384 + 0.02
0.461 0.017
0.011
0.015 + 0.0
Swiss roll high
1.0
Euclidean
0.349 0.004
0.409 0.003
0.0
0.018 = 0.0
Swiss roll very high
Diffusion Map
0.977 +: 0.005
0.984 0.004
0.018
0.026 + 0.0
Swiss roll very high
0.1
Heat-Geo
0.992 + 0.002
0.996 0.001
0.002
0.003 1 0.0
Swiss roll very high
0.1
Hcat-PHATE
0.991 + 0.001
0.997 0.001
0.079 0.003
0.101 = 0.003
Swiss roll very high
0.1
PHATE
0.631 + 0.01
0.594 = 0.011
0.023 1.0
0.028 - 0.001
Swiss roll very high
0.1
Rand-Goo
0.957 +: 0.002
0.994 0.001
0.009 1.0
0.012 0.0
Swiss roll very high
0.1
Shortest Path
0.999 1:: 0.0
0.999 0.0
0.006 0.0
0.007 + 0.0
Swiss roll very high
0.1
Euclidean
0.378 +: 0.002
0.121 0.002
0.0
0.018 - 0.0
Swiss roll very high
0.5
Diffusion Map
0.978 + 0.002
0,984 0.001
0.018
0.026 + 0.0 0.01 + 0.0
Swiss roll very high
0.5
Heat-Geo
0.997 + 0.0
0.998 0.0
0.00
Swiss roll very high
0.5
Heat-PHATE
0.996 +: 0.001
0.997 0.0
0.016
0.02 + 0.001
Swiss roll very high
0.5
PHATE
0.815 +: 0.002
0.823 0.004
0.032
0,049 + 0.002
Swiss roll very high
0.5
Rand- Geo
0.986 + 0.0
0.996 0.0
0.008 0.0
0.009 + 0.0
Swiss roll very high
Shortest Path
1.08:00.0 0.998
0.998 0.0
0.019 0.001
0.027 + 0.001
Swiss roll very high
0.5
Euclidcan
0.361 + 0.002
0.42 0.002
0.013 0.0
0.018 + 0.0
Swiss roll very high
1.0
Diffusion Mup
0.324 +: 0.061
0.399 0.033
0.018 + 0.0
0.026 + 0.0
Swiss roll very high
1.0
Heat-Geo-
0.466 +: 0.007
0.506 = 0.006
0.011 + 0.0
0.013 + 0.0
Swiss roll very high
1.0
Heat-PHATE
0.309 + 0.011
0.43 0.019
0.011
0.014 0.0
Swiss roll very high
PHATE
0.377 + 0.011
0.425 0.009
0.0
0.062 = 0.004
Swiss roll very high
1.0
Rand-Geo
0.398 + 0.009
0.516 0.008
0.01 0.0
0.012 0.0
Swiss roll very high
1.0
Shortest Path
0.367 + 0.018
0.413 0.016
0.0 0.0
0.015 0.0
Swiss roll very high
1.0
Euclidean
0.336 1 0.002
0.402 0.002
0.012 0.0
0.018 - 0.0
Table 6: Comparison of the estimated distance matrices with the ground truth geodesic distance
matrices on the Swiss roll dataset. Best models on average are bolded (not necessarily significant).
23

---

data
Noise level
Method
PearsonR
SpearmanR
Norm Fro N2
Norm inf N2
Tree
1.0
Diffusion Map
0.748 + 0.125
0.733 + 0.111
0.113 + 0.012
0.161 + 0.019
Tree
1.0
Heat-Geo
0.976 + 0.019
0.977 + 0.02
0.092 = 0.011
0.135 - 0.018
Tree
1.0
Heat-PHATE
0.918 1 0.032
0.885 1: 0.04
0.03 1 0.005
0.044 1: 0.007
Tree
1.0
PHATE
0.671 = 0.021
0.398 + 0.052
0.051 # 0.008
0.084 = 0.017
Tree
1.0
Rand-Geo
0.926 + 0.011
0,966 +: 0.019
0.076 + 0.01
0.117 + 0.018
Tree
1.0
Shortest Path
0.965 + 0.026
0.963 + 0.027
0.039 + 0.008
0.06 + 0.008
Tree
1.0
Euclidean
0.508 = 0.039
0.483 = 0.052
0.092 = 0.011
0.138 = 0.018
Treee
5.0
Diffusion Map
0.656 +: 0.054
0,653 +: 0.057
0.113 + 0.012
0.161 +: 0.019
Tree
5.0
Heat-Gen
0.822 +: 0.008
0.807 + 0.016
0.1 0.012 0.012
0.146 0.019
Tree
Heat-PHATE
0.765 0.025
0.751 = 0.023
0.043 0.006
0.08
Tree
PHATE
0.766 1: 0.023
0.743 1: 0.028
0.055 0.007
0.09
Tree
5.0
Rand-Geo
0.806 + 0.014
0.795 0.018
0.094 0.011
0.139 0.018
Tree
5.0
Shortest Path
0.78 = 0.009
0.757 0.019
0.075 0.009
0.117
Tree
5.0
Buclidean
0.735 1 0.014
0.704 1. 0.033
0.096 1 0.011
0.141 0.017
Tree
10.0
Diffusion Map
0.538 + 0.05
0,471 +: 0.089
0,113 + 0.012
0,161 + 0.019
Tree
10.0
Heat-Geo
0.62 0.025
+
0.I 0.012
0.146 0.019
Tree
10.0
Heat--HATE
0.63 1 0.018
0.588 0.031
0.046 1 0.005
0.083 1 0.012
Tree
10.0
PHATE
0.623 + 0.016
0.583 + 0.029
0.07 = 0.01
0.017
Tree
10.0
Rand-Geo
0.578 = 0.043
0.558 + 0.053
0.095 + 0.011
0.14
Tree
10.0
Shortest Path
0.539 +: 0.041
0.513 + 0.055
0.072 + 0.01
0.118 +: 0.017
Tree
10.0
Euclidean
0.508 = 0.039
0.483 + 0.052
0.092 + 0.011
0.138 + 0.018
Tree high
1.0
Diffusion Map
0.754 + 0.049
0.741 + 0.057
0.267 + 0.021
0.369 + 0.026
Tree high
1.0
Heat-Geo
0.996 = 0.001
0.999 + 0.001
0.242 = 0.02
0.338 = 0.026
Tree high
1.0
Heat-PHATE
0.927 1 0.011
0.875 1 0.032
0.062 1 0.003
0.084 0.006
Tree high
1.0
PHATE
0.528 + 0.085
0.141 0.061
0.209 0.023
0.027
Tree high
1.0
Rand-Geo
0.85 0.85+00.00 0.014
0.944 0.011
0.227 + 0.02
0.323 + 0.025
Tree high
1.0
Shortest Path
0.998 0.001
+.000
+
0.018 0.005
Tree high
Euclidean
0.928 1 0.018
0.928 1 0.024
0.24 0.02
0.334 1 0.026
Tree high
5.0
Diffusion Map
0.706 = 0.124
0.705 + 0.113
0.267 + 0.021
0.369 + 0.026
Tree high
5.0
Heat-Geo
0.97 : 0.01
0.975 0.009
0.253 1: 0.021
0.353 1: 0.026
Tree high
Hear-PHATE
0.932 = 0.022
0.919 0.03
0.072
0.112 0.008
Tree high
PHATE
0.913 1: 0.014
0.872 0.034
0.10
0,278 1: 0.025
Tree high
5.0
Rand-Geo
0.968
0.971
0.245 +-255+0.01
0.342 0.024
Tree high
5.0
Shortest Path
0.952 = 0.016
0.95
0.137 0.017
0.209 = 0.024
Tree high
5.0
Buclidean
0.882 1 0.028
0.873 1 0.032
0.237 0.02
0.333 1 0.025
Tree high
10.0
Diffusion Map
0.598 0.117
0.613 + 1.13
0.267 + 0.021
0.369 + 0.026
Tree high
10.0
Heat-Geo
0.861 0.039
0.87
0.254 0.021
0.353 + 0.026
Tree high
10.0
Heat--HHATE
0.844 1 0.05
0.838 0.055
0.168 0.015
0.27 1 0.025
Tree high
10.0
PHATE
0.837 + 0.052
0.838 + 0.049
0.204 0.018
0.301 + 0.024
Tree high
10.0
Rand-Geo
0.845 = 0.041
0.86 = 0.038
0.248 + 0.02
0.346 - 0.025
Tree high
10.0
Shortest Path
0.779 1 0.051
0.777 1 0.054
0.159 1 0.018
0.257 1 0.026
Tree high
10.0
Euclidean
0.709 + 0.054
0.699 + 0.059
0.229 + 0.02
0.327 + 0.026
Table 7: Comparison of the estimated distance matrices with the ground truth geodesic distance matrices on the Tree roll dataset. Best models on average are bolded (not necessarily significant).
C.7 Time Complexity
present the average computing time for creating embeddings and corresponding distance matrix for the different methods. All methods are applied on the Swiss roll dataset in three dimension with 2000 samples. We present empirical averages and standard deviations over ten repetitions. The experiments were run on a Apple M2 Pro chip with 16G RAM.
In Table 12
24

---

Nois: level
Method
PearsooR
SpearmanR
Norm Fro N2
Norm N2
Diffusion Map
0.974
0.985
Swiss
Hear-Gion
0 0.003
0.99 0.002
0
0.026 0.0
Swiss nill
0.1
Heat-PIAATE
0.00 1: 0.002
0.907
0.001
0.018 1 0.0
0.026 0.0
Swiss mlll
0.1
PHATE
0.677 0.02
0.607 0.014
0.015+ + 0.0
0.026 + 0.0
Swiss roll
0.1
Rand-Geo
0.917 1 0.003
0.915 1 0.002
0.018 1 0.0
n.1226 1. 0.0
Swiss roll
0.1
Shanest Puth
1.000.0 0.0
1.000.0
0.018 + 0.0
0.026 = 0.0
Swiss
TSNE
0.905
0.897
0.006 2.0.
0.088 0.0
Swiss
UMAP
0.013
0.70
0.00 0.018
0.016 0.001
Swiss
Euclidean
0.003
0.421 0.003
Swiss
Diffusion Map
0.992 0.003
0.957 0.002
0.015
0.026 0.0
Swiss ril
0.5
Hear-Geo
0.997
0.1
0.39% + 0.001
0.018 + 11.
0.026 0.0
Swiss nill
0.5
Hen-PEIATE
0.933 2: 0.001
0.907 : 0.0
0.018 :: 0.0
0.026 1: 0.0
Swiss roll
0.5
PHATE
0.696 2: f10111
0.711 = 0.008
0,018 4: 11.
IJ126 2 0.0
Swiss rull
Gend-Geo
0.952 255220.000
0.932 0.002
0.018 1
0.026 1 0.0
Shantest
0.0
99999
0.018 +
0.0
Suiss
TSNE
0.899
0.802
Suiss
UMAP
0838+0.010
0.810 0.017
0.012
0.016 0.001
Swiss roll
0.5
Euclidean
0.581 0.004
0.421 0.004
0.018 1 0.0
0.026 0.0
Swiss mil
1.0
Diffusion Map
0.476 0.226
0.478 0.138
0.018 + 0.0
0.026 + 0.0
Swist roll
1.0
Hear-Geo
(.072 11:221
0J676 0.193
0.018 + 1LO
DJI22 = 0.0
Swiss nill
1.0
Heat-PRAATE
0.674 0:::::0..000 0.160
0,684 0.154
0.018 + 0.0
0.026 + 0.0
Swiss
PHATE
0.287 0.00
0.349 0.028
0.018 1:
IL1126 0.0
Rand-Gon
0.029
0.022
0.018
Swiss
Shortest
0.467 0:17
0.511 0.163
0 +
0.0
Swiss
TSNE
4.118+0.188
0.724 0,151
O 000 0.002
0.044 003
Swiss mill
1.0
UMAP
0.727 + 0.181
0.713 + 0.167
0.012+ 0.001
0 0.001
Swiss nill
1.0
Euchideam
0.571 1 0.006
0.414 1 0.005
0.018 1 0.0
0.026 0.0
Swiss rell
5.0
Diffusion Map
0.157 + 0.021
0.173 + 0.015
0.018 + 0.0
0.026 + 0.0
Swiss nill
5.0
Hea-PPAATE
0.235 0.014
0 239 1 0.013
0.018 1: 0.0
0.026 1. 0.0
PHATE
0.201
0.014
0.237
Rand-Geo
0.201
0.238
0.018 1.
Swiss
Shortest
0.233 0.01
0.018
0.026
Swiss roll
5.0
TSNE
02200.111
0,233 = 0.01
0.012 0.0
0.018 0.00000000 0U0
Swiss mill
5.0
UMAP
0:205 + 0.013
0. 0:239
0.012
0.015 0.0
0.022 0.0
Swiss roll high
0.1
Diffusion Map
(138 2: 0.003
0.986
0
0.018 =: 1L.
ILI22 d: 0.0
Swiss roll high
0.1
Heat-Gee
0.906 + 0.002
0 0.001
0.018 + 0.0
0.026 + 0.0
Swiss high
Heat-PHATE
0,991
0.002
0.018
high
PHATE
0.678
0.608
0.015
Swiss bigh
Rand-Geo
0.917 0.003
0.915
0.018 1:
Swiss roll high
0.1
Shortest Puth
1.0 0.0
1.0 1 0.0
0.018 0.0
0.026 0.0
Swiss roll high
0.1
TSNE
0.903 0.000
0.89 = 0.003
0.00G 0.0
D 0.0 = 0.0
Swiss nill high
011
UMAP
0.806 +: 0.014
0.794 + 0.01
0.011 + 0.0
0.016 +: 0.001
Suiss roll high
0.1
Fuclidean
0.394 + 0.002
0.424 0.002
0.018 + 0.0
0.026 0.0
Swiss roll high
Diffusion Map
0.002
0.985 0,002
0,018 1:
IJI26 0.0
Swiss high
Heat-Gea
0.938 + 0.0
0.907
0.018 +
0.026
Swiss high
Hear-PHATE
0.95= 0.0
0.997
0.018
Swiss high
PHATE
0.754
0.01
0.756 0.006
0.018 +
0.026 0.0
Swiss roll high
0.5
Rand Geo
0.945 0.001
0,945 0.002
0.018 1: 0.0
0.026 1 0.0
Swiss zoll high
0.5
Shortest Path
0.999 0.0
0 998 1 0.0
0.018 1. 0.0
0.026 1 0.0
Swist roll ligh
0.5
TSNE
(L9G5 11006
0.899 0.003
0.00G 0.0
D 0ON 0.0
Swiss rill high
0.5
UMAP
0.876 3776+-0.011 0.017
0.86 - 0.024
0.012 0.0
0.017 + 0.001
high
Fuelidean
0.000
bigh
Diffusion Map
0.155
0.526
0.018
Swiss high
Hear-Gion
0.643 0.173
0,633
0.018
0.0
Swiss nill high
1.0
Hea-PPHATE
0.609 0.17
0.611 0.121
0.018 0.0
0.026 0.0
Swiss mill high
1.0
PHATE
0.271 0.025
0.343 0.011
0.018 0.0
0.026 0.0
Swiss zull high
1.0
Ran
0.41 0.028
0.446 1: 0.08
0.018 1 0.0
IJ1266 1 0.0
Swiss roll ligh
1.0
Shartest Path
0.343
0.4 = 0.007
0 + 0.0
B.I26 = 0.0
Swiss high
TSNE
0.737 0.1124
0.725
0.008
0.0U3
Swiss high
UMAP
0.803
860 0.083
+ 0.001
0.02 0.001
Swiss bigh
Euclidean
0.003
0.414 0.002
0.018
Swiss high
Diffusion Map
0.161+0.011
0.174 0.009
0.015
0.026 0.0
Swiss roll high
5.0
Hear-PHATE
0.202 0.01
0.236 + 0.009
0.018 (1.
0. 0.0
Swiss noll high
5.0
PHETTE
0.201
0.01
0.251 :: 0.008
0.018 0.0
0.026 :: 0.0
Swiss zoll high
5.0
Ran
0.192 1L009
0.228 0.008
0,018 2: 11.
I.1266 2: 0.0
Swiss zull high
5.0
Shortest Path
0,187 1.001
0.221 0.009
0.018 1. 0.0
0.026 1 0.0
Swiss ligh
TSNE
0.182 0.011
0.213 0.01
0.013=0.0
0.019
high
UMAP
0.195
0.016::0.0
0.001
Swiss roll very high
0.1
DifTusion Map
0.977 0.0000000
0.984 0.004
0.018 4: 0.0
0.026 0.0
Swiss roll very high
0.1
Hear-Geo
0.996 1 0.000
0.997 0.001
0.018 1: 010
0.026 1: 0.0
Swiss mill very high
0.1
Hea--PHATE:
0.901 + 0.001
0.907
0.001
0.018 0.0
0.026 + 0.0
Swiss roll very high
0.1
PHATE
(1883 2 I1023
0.701
0.016
0.018 = 1.0
DJ126 = 0.0
Swiss nill very high
0.1
Rand-Gioo
0 0.002
0.917 0.002
0.018 + 0.0
0.026 0.0
Swiss very high
Shortest Path
0.0
0.99942
0,018 1.
1L1126 0.0
Swiss very high
0,897 0.004
0.011 = 1.0.0
Swiss very high
UMAP
0U785 0.114
0.781 0.117
0.016 0.001
Swiss mill very high
0.1
Fuclidlea
0.334 + 0.002
0.424
0.002
0.018 +
0.026 0.0
Swiss roll very high
0.5
Didfusion Map
0.978 2: 0.002
0.984 = 0.001
0.018 (LO
IL =: 0.0
Swiss zoll very high
0.5
Hear-Geo
0.997
0.998 100
0,018 1 0.0
IJR266 1 0.0
Suiss rlll very high
0.5
Hear-PHATE
0.996 0.090=0..00 (1.01
0.907 0.0
0.018 + 0.0
D.826 = 0.0
Swiss nill very high
PHETTE
0.827 1.7.1.00000 0.005
0.815 0.002
0.018 1:
0.026 0.0
very high
Rand-fino
0.944+00.00
0.944
very high
Shortest
0.998 0.0
0.997
0.018
Swiss very high
0.177
0.97 0.007
O. 0.000 0.0
0.008 0.001
Swiss roll very high
0.5
UMAP
0.928 1.928=0000 0.011
0.925 = 0.012
0.012 (LO
0 = 0.001
Swiss mill very high
0.5
Fuclidian
0.579 + 0.002
0.42 - 0.002
0.015 + 0.0
0.026 + 0.0
Swiss zell very high
1.0
Diffusion Map
(L324 (L06]
0,399 = 0.033
0.018 ILO
ILR26 A: 0.D
Suiss mill very high
1.0
Heat-Gien
0.364 0.008
0.425 + 0.015
0.015 + 0.0
0.026 + 0.0
very high
Heat-PHATE:
0.352
0.018
0.0
very high
PHATTE
0.326 0.000
0.388
0.015 0.0
Swiss very high
Rand-Gor
0.357 0.007
0.604
0.018
Swiss roll very high
1.0
Shortest Pash
0.335 0.001
0.39 L 0.011
0.018 1 0.0
0.026 0.0
Swiss roll very high
1.0
TSNE
0.515
0.522 = 0.01
0.012=00.0
0.016 0.0
Swiss nill very high
1.0
UMAP
0.705 : 0.059
0 0.787 : 0.058
0.015
0.021 : 0.0
Suiss woll very high
1.0
Fuclideean
0.386 + 0.002
0.41 - 0.001
0.015 + 0.0
0.026 0.0
Table 8: Comparison of the estimated distance Ahatiiess with the ground truth geodesic distance
matrices on the Swiss roll dataset, using two-dimensional embedding. Best models on average are bolded (not necessarily significant).

---

data
Noise level
Method
PearsonR
SpearmanR
Norm Fro N2
Norm N2
Tree
0.1
Diffusion Map
0.748 + 0.125
0.733 + 0.111
0.113 + 0.012
0.161 + 0.019
Tree
0.1
Hcat-Geo
0.943 + 0.037
0.94 + 0.037
0.113 = 0.012
0.161 # 0.019
Tree
0.1
Heat-PHATE
0.872 1 0.04
0.83 1 0.061
0.113 1 0.012
0.161 1 0.019
Tree
0.1
PHATE
0.564 + 0.039
0.469 + 0.052
0.113 0.011
0.161 + 0.018
Trree
0.1
Rand-Gco
0.937 0.868 = 0.017
0
0.113
0.161 0.019
Shortest Path
0.931 0.041
0.012
0.161 0.019
Tree
TSNE
0.847 + 0.034
0.045
0.123 +
Tree
0.1
UMAP
0.692 + 0.058
0.671 +0047
0.107+0.00
0.153+0.015
Tree
0.1
Euclidean
0.809 = 0.017
0.778 = 0.024
0.113
0.011
Tree
Diffusion Map
0.656 + 0.054
0.653 + 0.057
0.113 0.012
0.161 + 0.019
Tree
0.5
Heat-Geo
0.806 0.019
0.787 0.009
0.113 0.012
0.161 0.019
Tree
0.5
Heat-PHATE
0.746 1: 0.024
0.744
0.012
0.161 0.019
Tree
0.5
PHATE
0.766 = 0.023
0.746
0.011
0.018
Tree
0.5
Rand-Geo
0.721 +: 0.024
0.024
0.113
0.019
Tree
0.5
Shortest Path
0.765 + 0.01
0.738 0.011
0.113
0.161 0.019
Tree
0.5
TSNE
0.795 = 0.046
0.766 0.055
0.083
0.128
Tree
UMAP
0.783 = 0.06
0.757 0.054
0.11
0.157
Tree
Euclidean
0.704 0.02
0.672 + 0.038
0,113
0.161 0.019
Tree
1.0
Diffusion Map
0.338 = 0.05
0.471 # 0.089
0.113
0.161
Tree
1.0
Heat-GGeo
0.613 1 0.025
0.58 0.0066
0,113 0.012
0.161
Tree
1.0
Hear-PHATE
4.14-00002
0.571+0.044
0.012
+
Tree
PHATE
0.615
0.036
0.11
0.015
Rand-Geo
0.487 0.064
0.0071
0.113
0.019
Tree
Shortest Path
0.542 0.047
0.06
0
0.019
Tree
1.0
TSNE
0.583 +: 0.042
0.553 0.045
0.086 1 0.011
0.135 - 0.017
Tree
1.0
UMAP
0.595 = 0.032
0.562 = 0.036
0.111 = 0.011
0.158 0.019
Tree
1.0
Euclidean
0.502 1 0.051
0.479 1 0.064
0.113 1 0.012
0.161 1 0.019
Tree high
0.1
Diffusion Map
0.754 0.049
0.741 +: 0.057
0.267 1 0.021
0.369 1 0.026
Tree high
0.1
Heat-Geo
0.956 + 0.014
0.957 + 0.015
0.267 + 0.021
0.369 + 0.026
Tree high
0.1
Heat-PHATE
0.831 = 0.082
0.764 0.155
0.267 0.021
0.369
Tree high
0.1
PHATE
0.484 1 0.036
0.410.022
0.267 0.02
0.369 1 0.025
Tree high
0.1
Rand-Geo
0.817 + 0.013
4774+00022
0.267 + 0.021
0.369 + 0.026
Tree high
0.1
Shortest Path
0.958 + 0.014
0.956 0.017
0.267 + 0.021
0.369 + 0.026
Tree high
0.1
TSNE
0.89 0.039
0.866 0.044
0.233 1 0.021
0.327 0.026
Tree high
UMAP
0.8
0.764
0,259 + 0.021
0.36 0.028
Tree high
Fuclidean
0.878 0.042
0.859 + 0.051
0.267 + 0.021
0.369 0.026
Tree high
0.5
Diffusion Map
0.706 1 0.124
0.705 1 0.113
0.267 1 0.021
0.369 1 0.026
Tree high
Heat-Geo
0.932 0.022
0.928 0.023
0.267 + 0.021
0.369 + 0.026
Tree high
Heat-PHATE
0.923 0.023
0.921 + 0.022
0.267 + 0.021
0.369 - 0.026
high
PHATE
0.844 0.048
0.79 0.07
0.267 1
0.369 1 0.025
Tree high
Rand-Gco
0.875 0.042
0.855 0.048
0.267 = 0.021
0.369 = 0.026
Tree high
0.5
Shortest Path
0.917 0.025
0.91 0.03
0.267 1: 0.021
0.369 1 0.026
Tree high
0.5
TSNE
0.922 + 0.035
0.91 + 0.045
0.237 + 0.021
0.334 0.027
Tree high
0.5
UMAP
0.823 = 0.054
0.803 = 0.041
0.261 = 0.021
0.361 = 0.026
Tree high
0.5
Euclidean
0.819 + 0.048
0.799 + 0.053
0.267 + 0.021
0.369 + 0.026
Tree high
1.0
Diffusion Map
0.598 + 0.1
0.613 + 1.33
0.267 + 0.021
0.369 + 0.026
Tree high
1.0
Heat-Geo
0.794 = 0.066
0.805 + 0.049
0.267 + 0.021
0.369 + 0.026
Tree high
1.0
Heat-PHATE
0.826 1 0.064
0.823 1 0.067
0.267 1 0.021
0.369 1 0.026
Tree high
1.0
PHATE
0.827 + 0.059
0.82 1 0.062
0.267 + 0.02
0.369 + 0.025
Tree high
1.0
Rand-Geo
0.71 + 0.043
0.686 + 0.045
0.267 + 0.021
0.369 + 0.126
Tree high
1.0
Shortest Path
0.771 0.064
0.753 1 0.07
0.267 1 0.021
0.369 1 0.026
Tree high
1.0
TSNE
0.84 0.066
0.821 + 0.074
0.238 0.02
0.335 - 0.026
Tree high
UMAP
0.853 0.051
0.839 0.057
0.264 1 0.021
0.365 1 0.026
Tree high
1.0
Euclidean
0.683 0.067
0.665 0.07
0.267 0.021
0.369 0.026
Table 9: Comparison of the estimated distance matrices with the ground truth geodesic distance matrices on the Tree dataset, using two-dimensional embedding. Best models on average are bolded (not necessarily significant).
26

---

data
Noise level
Method
Homogeneity
Adjusted Rand Score
Adjusted Mutual Info Score
Swiss roll
0.1
Heat-Geo
0.82 + 0.008
0.668 = 0.034
0.74 0.018
Swiss roll
Phate
0.731 0.035
0.546 = 0.044
0.652 = 0.046
Swiss roll
0.1
TSNE
0.748 + 0.067
0.537 + 0.1
0.668 + 0.068
Swiss roll
0.1
UMAP
0.81 + 0.036
0.611 + 0.039
0.726 + 0.045
Swiss roll
0.5
Heat-Geo
0.813 0.026
0.656 = 0.049
0.733 = 0.022
Swiss roll
0.5
Phate
0.735 + 0.048
0.543 + 0.064
0.656 + 0.053
Swiss roll
0.5
TSNE
0.764 0.07
0.564 + 0.097
0.684 +: 0.065
Swiss roll
UMAP
0.826 0.019
0.664 0.073
0.744 0.032
Swiss roll
1.0
Heat-Gen
0.722 0.051
0.548 + 0.091
0.652 + 0.056
Swiss
Phate
0.482 1: 0.014
0.317
0.428 1: 0.021
Swiss roll
1.0
TSNE
0.757 0.037
0.562 0.058
0.679 0.042
Swiss roll
1.0
UMAP
0.041
0.51 = 0.077
0.65 + 0.05
Swiss roll high
0.1
Heeat-Geu
0.82 1: 0.015
0.666 1 0.033
0.739 1 0.019
Swiss roll high
0.1
Phate
0.705
0.518 + 0.048
0.628 + 0.04
Swiss roll high
0.1
TSNE
0.558 = 0.115
0.677 + 0.08
Swiss roll high
0.1
UMAP
0.796 0.03
0.624 + 0.048
0.714 = 0.037
Swiss roll high
0.5
Heat-Geo
0.805 + 0 0.021
0.655 = 0.047
0.725 1 0.035
Swiss roll high
Phate
0.745
0.562 = 0.061
0.664 = 0.047
Swiss roll high
TSNE
0.747
0.538 0.11
0.668 + 0.075
Swiss roll high
0.5
UMAP
0.787+
0.573 + 0.067
0.703 + 0.032
Swiss roll high
Heat-Gco
0.045
0.534 0.057
0.644 0.032
Swiss roll high
1.0
Phate
0.552 0.047
0.386 + 0.056
0.496 + 0.04
Swiss roll high
1.0
TSNE
0.754 4
0.548 + 0.068
0.675 + 0.036
Swiss roll high
1.0
UMAP
0.76 1: 0.041
0.56 1 0.077
0.68 1 0.05
Swiss roll very high
0.1
Heat-Geo
0.818 + 0.033
0.668 + 0.074
0.738 + 0.039
Swiss roll very high
0.1
Phate
0.688 1: 0.043
0.497 1 0.053
0.614 1: 0.053
Swiss roll very high
TSNE
0.007
0.544 1 0.101
0.662 1 0.075
Swiss roll very high
0.1
UMAP
0.016+0.042
0.65 0.069
0.733 + 0.054
Swiss very high
Heeat-Gev
0.73 0.0045
0.605 0.093
0.701 1: 0.028
Swiss roll very high
0.5
Phate
0.758.10.034
0.55 0.037
0.676 0.014
Swiss roll very high
0.5
TSNE
0.05
0.557 : 0.093
0.708 + 0.031
Swiss roll very high
0.5
UMAP
0.789 1: 0.052
0.574 1 0.101
0.707 1: 0.061
Swiss roll very high
1.0
Heat-Geo
0.592
0.427 + 0.063
0.545 - 0.031
Swiss roll very high
1.0
Phate
0.377 + 0.046
0.486 + 0.045
Swiss roll very high
TSNE
0.738 0.019
0.551 0.039
0.662 0.025
Swiss roll very high
UMAP
0.736 0.057
0.542 1 0.102
0.66
Tree
0.1
Heat-Geo
0.784 0.051
0.734 - 0.07
0.786 0.051
Tree
Phate
0.042
0.409 0.064
0.555 + 0.042
Tree
TSNE
0.706
0.61 0.075
0.712 1 0.055
Tree
0.1
UMAP
1678
0.584 + 0.12
0.681 = 0.086
Tree
0.5
Heat-Gee
0.
0.411 + 0.154
0.577 + 0.094
Tree
0.5
Phate
0.529 0.
0.404 + 0.151
0.555 1: 0.095
Tree
0.5
TSNE
0.647 + 0.049
0.591 + 0.065
0.65 0.048
Tree
0.5
UMAP
0.645
0.565 + 0.058
0.652 0.05
Tree
1.0
Heat-Geo
0.398 0.07
0.3 0.077
0.42 1 0.07
Tree
1.0
Phate
0.418 0.08
0.337 + 0.093
0.43 + 0.075
Tree
TSNE
0.405 0.077
0.378 0.074
0.405 + 0.077
Tree
1.0
UMAP
0.432 0.086
0.395 1 0.098
0.432 0.085
Table 10: Clustering results on swiss roll (with distribution) and tree. Best models on average are bolded (not necessarily significant).
27

---

Embeddings of HeatGeo for different hyperparameters
Harnack Regularization
au2
WMM
-C-1
-AU2
-10S
Harnack 0
Hamack - U.
Hamack
Harnack 1.5
Hamack
Diffusion Time
-000
-192
0M0
-105
-081
ME
-422
-201
-000
-291
Que
8.31
-101
6ill
Ciffusion meme
Difluion Time
DiRusio Time
10
Difindon Time 20
Diftason Time / wuto
Number of neighbours
E3
0/0
0:11
0/1
000
3801
0.0
3101
0213
000
8.0 000
6:11
20
Order of Heat Kernel Approximation (Euler)
0.82
1.32
one
100
0111
6.30
01Z
200
DRO
8:01 OLL
8:11
000
Order Heat Kernel Approximation (Chebyshev)
out
1022
-000
(U0
Figure 8: Embeddings of Heat Geodesic Embedding for different choices of hyperparameters on the dataset. We evaluate the impact of the Harnack regularization, the diffusion time, the number of neighbours in the kNN, and the order of the approximation for Euler and Checbyshev approximations.
28

---

Pearson-p in function of on Swiss Roll data
1.0
Configuration: Dataset:Swiss Roll Noise:0.5
0.8
Pearson-p
0.6
Heat Geo
0.4
Heat
Heat PHATE
Heat PHATE
Rand Geo
0.2
Phato
Phate
Diffusion Map
20
40
60
80
100
Figure 9: Impact of diffusion time on the Pearson correlation between the estimated distance matrix and ground truth distance matrix for different methods on the Swiss roll dataset.
Pearson-p in function of Order of the heat kernel approximation on Swiss Roll data
1.0
Hest Geo
PHATE
Configuration: Dataset Swiss Roll
0.8
Noisc:0 5
Filter Methoc/Cheebssher Approximation
Pearson-p
0.6
0.4
0.2
20
Order the heat kernel approximation
Figure 10: Impact of Checbyshev approximation order the Pearson correlation between the estimated distance matrix and ground truth distance matrix for different methods on the Swiss roll dataset.
29

---

Pearson-p in function of k on Swiss Roll data
1.0
Configuration: Datasct:Swiss Roll No.se:0.5
0.9
Pearson-p
0.8
0.7
Heat Geo
Heat PHATE
Rand Geo
0.6
Phate
Diffusion Map
Shortest Path
10
15
20
25
30
Figure Impact of number of neighbours on the Pearson correlation between the estimated distance matrix and ground truth distance matrix for different methods on the Swiss roll dataset.
Pearson-p p in function of Harnack Regularization on Swiss Roll data
0.99
Configuration:
Dataset:Swiss Roll
0.98
Moise:0.5 Filter Method:Culer Approximation
Pearson-p
7:auto
0.97
0.96
0.95
0.00
0.25
0.50
0.75
1.00
L.25
1.50
Harnack Regularization
Figure 12: Impact of Harnack regularization the Pearson correlation between the estimated distance matrix and ground truth distance matrix for HeatGeo on the Swiss roll dataset.
30

---

Pearson p in function of Graph type on Swiss Roll 10D data
Heat Geo
0.900
0.898
0.896
Pearson p
0.894
0.892
0.890
0.888
alpha
knn
scanpy
umap
Graph type
Figure 13: Pearson correlation between estimated and ground truth distance matrices for the 10- dimensional Swiss roll dataset for various graph constructions. Standard deviations are computed the test folds.
Table 11:
Clustering quality metrics for different methods. We report the homogeneity and the adjusted mutual information (aMI), Best models on average are bolded (higher is better).
data
Noise level
Method
Homogeneity
Adjusted Rand Score
Adjusted Mutual Info Score
Mnist
0
Diff-Map
0.556 + 0.002
0.347 0.002
0.622 = 0.002
Mnist
Heat-Geo
0.785 + 0.0
0.695 + 0.0
0.829 + 0.001
Mnist
Phate
0.822 + 0.01
0.72
0.835 + 0.011
Mnist
TSNE
0.903 = 0.003
0.871 0.002
0.902 = 0.003
Mnist
UMAP
0.851 + 0.016
0.846 + 0.005
0.86 0.015
Coil
Diff-Map
0.21 0.036
0.041 + 0.015
0.142 + 0.024
Coil
Heat-Geo
0.849 + 0.016
0.67 0.029
0.806 + 0.022
Coil
0
Phate
0.804 + 0.017
+ 0.615 + 0.028
0.735 + 0.021
Coil
0
TSNE
0.907 + 0.014
0.79 0.03
0.88 1: 0.02
Coil
0
UMAP
0.871 + 0.009
0.826 + 0.012
Pbme
0
Diff-Map
0.026 + 0.001
0.011 + 0.0
0.038 # 0.001
Pbme
0
Heat-Geo
0.734 + 0.009
0.721 + 0.019
0.768 + 0.017
Pbmm
Phate
0.798 = 0.012
0.818 0.009
0.785 0.01
Phme
TSNE
0.605 + 0.019
0.437 + 0.032
0.544 + 0.022
Phme
UMAP
0.177 + 0.037
0.097 + 0.033
0.148 + 0.035
Table 12: Average computing time for creating embeddings and corresponding distance matrix for the different methods. All methods are applied on the Swiss roll dataset in three dimension with 2000 samples. We present empirical averages and standard deviations over ten repetitions.
Method
Time (s)
UMAP
4.21 0.68
I-SNE
80.15 62.44
Isomap
2.58 + 0.03
PHATE
3.30 0.05
Diff Map
8.62 1: 0.23
HeatGeo (Backward Euler)
2.60 + 0.05
HeatGeo (Chebyshev)
2.11 0.05
31


Reviews:
[{"review_id": "ZfBwcwoTEK", "reviewer": "Reviewer_7Dv3", "summary": "This article deals with representation learning and dimensionality reduction of high dimensional noisy datasets. The proposed method relies on heat diffusion. More precisely, assuming their data can be viewed as Riemannian manifolds and based on Varadhan’s formula, the authors define the notion of heat-geodesic dissimilarity. According to the authors, using the heat equation guarantees the conservation of geodetic distances between the raw data and its integrated counterpart. This dissimilarity is then used through multidimensional scaling (MDS) to define their so-called heat geodesic embedding.\n\nAfterward, the authors theoretically compare their method to the start-of-the-are ones. In particular, in the appendix, they attempt to create a generalized framework for dimensionality reduction methods. \n\nLast, they conduct numerical experiments to test their algorithm on different data types. In particular, they want to ensure that their method effectively preserves geodetic distances. They also seek to ensure that the inherent data structures are conserved: Preserving clusters or data temporality, for instance.", "questions": "Some comments:\n- Line 66: $n$ not defined.\n- Line 67: I think it should be $\\frac{n}2$ instead of $2n$.\n- Eq (2): $V(x,\\sqrt{t})$ not defined.\n- Paragraph starting line 79 (Graph construction and diffusion): $n$ is not defined, although we understand it is the sample size in context.\n- Line 88: Notation $[n]$ not defined.\n- Line 94: Add references when referring to literature.\n- Line 160: Should the $n$ not be a $d$ here?\n- Line 181: they use the a distance... -> the _or_ a.\n- Line 304: Strange line break.", "limitations": "Not concerned.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- Interesting idea, based on the solid formalist of Riemannian manifolds.\n- Efforts to unify different embedding approaches (especially in the appendix) with comparisons with existing methods.\n- Numerical experiments demonstrating the benefits of heat geodesic embedding for understanding data.\n", "weaknesses": "There are typos and undefined notations in several places (see below). This is a disservice to the paper.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their suggestions, and we agree that the quality of the writing could be improved. We did a thorough pass on the manuscript, and corrected some typos including the ones mentioned. We also added the missing definitions. We believe that these minor corrections help improve the overall clarity of the manuscript, and respectfully ask the reviewer to reconsider their impression of the paper and potentially improve the given score if the raised concerns have been allayed. We thank the reviewer again for their time and we are also happy to answer any further questions that arise!"}, {"review_id": "3cDJlEmQk6", "reviewer": "Reviewer_B5m7", "summary": "The paper introduces a dimensionality reduction technique based on heat kernels and called HeatGeo. The main idea is to learn a mapping that projects points to a low-dimensional space (ideally 2 or 3-dimensional for visualization) so that relative distances between pairs of points satisfy the same order as distances on some (complete Riemannian) manifold that uses some distance based on the heat kernel. ", "questions": "I have no questions.", "limitations": "N/A", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Originality: The idea of using heat kernels for dimensionality reduction is not new (e.g. Belkin & Niyogi [1] that proposed Laplacian eigenmaps). This paper proposes a generalization of previous techniques by exploiting Varadhan's formula (Eq. (1)) and by bounding the calculated distance on a low-dimensional complete Riemannian manifold. They show that the ordering of distances in the high-dimensional space tends to be preserved on the low-dimensional manifold.\n\nQuality and Clarity: The paper is well-written and mentions different approximation techniques that can be used to approximate the computation of the kernel matrix since the algorithmic complexity is cubic in the number of points. Overall, the paper shows that the proposed approach tends to preserve clusters and the ordering of distances in the low-dimensional space better than standard approaches such as t-SNE and UMAP.\n\nSignificance: The significance of the approach is difficult to quantify. As most dimensionality reduction techniques that are often used for qualitative evaluation in the machine learning literature, the proposed approach can improve the understanding of a data set with a visualization that preserves the original general structure.", "weaknesses": "I do not see a particular weakness in the formulation of the approach which has a nice theoretical justification. \n\nThe main weakness of the paper is that it mainly focuses on dimensionality reduction and visualization which already have widely used techniques (e.g. t-SNE and UMAP) that are used only for qualitative evaluation in practice.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their encouraging review. Below we respond to the main weakness raised by the reviewer.\n\n>The main weakness of the paper is that it mainly focuses on dimensionality reduction and visualization which already have widely used techniques (e.g. t-SNE and UMAP) that are used only for qualitative evaluation in practice.\n\n\nWe agree with the reviewer that several methods have been proposed for data representation learning. A major motivation for our work is to unify these different methods through the lens of differential geometry (with Varadhan formula and Harnack inequality). We believe this helps understand why and when these methods work and how their algorithmic difference can impact the resulting embeddings. For instance, neither t-SNE nor UMAP preserve distances, and this has an impact on the visualization as opposed to methods that use diffusion or other discretizations of heat flow on the manifold which is explicitly connected to heat diffusion. \n\nWhile these methods can be assessed qualitatively, by inspecting the low dimensional embeddings, we aimed at also providing a rigorous quantitative evaluation. We did so by comparing the estimated distance matrices against ground truth, and by assessing whether the embeddings conserved the underlying structure of the data (through cluster evaluation and EMD interpolation). These results showed that our method, supported by its strong theoretical foundations, along with closely related approaches (e.g. PHATE), produced embeddings that were more faithful to the inherent structure of the data.\n"}, {"review_id": "21BuyFFL7L", "reviewer": "Reviewer_QP8B", "summary": "The authors develop specific heat diffusion kernels for manifold embedding.\nThey define a heat kernel based on the Varadha formula.\nThe authors try to weave some connections with existing methods of DR such as PHATE, t-SNE, etc.\nThey run some experiments", "questions": "What is V(x,sqrt(t)) ? please define and/or give an intuitive interpretation. The section about Varadhan's formula seems to end prematurely. Probably Vradhan could gain to be detailed more while the subsequent paragraphs in Preliminaries could be shortened or dropped (except for the quick reminder about stress-based MDS).\nSection 3 about related work is going too fast and lacks details (reminder of P for PHATE? drop SNE that is quite well known).\nIn section 4, if distance order preserving matters, why not compare to non-metric MDS too?\nCan you clarify Section 4? Definition 4.1 is cryptic and unintuitive; what is the meaning of sigma?\nIn proposition 4.1, what is |x-y| ? absolute value? distance like ||.||_2 elsewhere?\nWhy come back on MDS? Merge the presentation of MDS and its specific weighting?\nSection 5 is also wanting to embrace everything within limited space: alonger version in appendix and more details about the core contribution of the paper\nRef 33: Cédric Villani has a twin brother?", "limitations": "Yes. They have. No potential negative societal impact here.", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The paper is well written except for the organization, clarity, and self-containedness.\nThe presented developments about the geodesic distance and the heat diffusion kernel are of interest for the manifold learning community.\nThe authors quite honestly acknowledge some of the limitations of their work.", "weaknesses": "The writing style is too dense, poorly structured, with a succession of (too) short sections, referring to each other too many times, and often missing key information.\nIsomap is not cited and not compared to (in Fig. 1 for instance; later it seems that Isomap is actually used, without saying the name and still without a bibliograpcal reference, or is 'shortest path' with stress-based MDS instead of classical MDS? This can be referenced too, most probably).\nThe mathematical notation is misleading and unusual for manifold learning (x and y in the same space! use x_i and x_j ?).\nThe experiments have unclear/contradictory objectives: to detect manifold structures and clusters? The COIL-20 or 100 data set should then be considered. The results in Fig. 1 are not convincing at all: Isomap does better. In experiments related to clustering, the reader might always question the purposefulness of assessing cluster separation in methods which are local and depending on graphs; a disconnected graph somehow automatically guarantees cluster separation. On the other hand, the arrangement of clusters with respect to each other might actually be more relevant. The clustering experiments on Swiss roll and tree are far-fetched (use COIL-20?), with completely artificial clusters. Over-expression of clusters or even spurious clusters have been reported with many 'local' methods of manifold learning. Assessing the preservation of distances or neighbourhoods might be less of a stretch between the focus on geodesic distances and the current interest for clustering with those methods in communities of practitioners.  \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewers for their constructive comments on our manuscript. Below we will address the main concerns of the reviewer.\n\n>The writing style is too dense, poorly structured, with a succession of (too) short sections, referring to each other too many times, and often missing key information.\n\nThanks for this observation, we will try to improve the connection between paragraphs to make reading smoother.  See also our comments in the section **On the reviewer’s questions** for more detailed modifications on the organization of the manuscript.\n\n**About Isomap**\n\nIn our submission, in particular in Tab.1, we referred to Isomap as the shortest-path method since we only used the distance matrix (i.e. only the shortest path). Indeed, Isomap is MDS done on shortest path distances between points on an affinity graph. However, we improved the clarity of the paper by linking to the shortest path distance to Isomap, and referencing Isomap explicitly. We emphasize that while Isomap indeed finds manifold distances in noiseless settings, the performance degrades in noisy data where the data graph is not reliable. For instance, in the embryoid body data, Isomap muddles many of the differentiation structures as also noted in the PHATE paper [Moon et al. 2019] (see Fig.1 in the global response's PDF). The same observation is made in our new quantitative results on Isomap (see tables in the global response's PDF)\n\n**About misleading notation**\n\n>The mathematical notation is misleading and unusual for manifold learning (x and y in the same space! use x_i and x_j ?). \n\nWe have uniformized the notation, we use $x_i,x_j$ for points in the ambient space and $y_i,y_j$ for points in the embedded space.\n\n**About the apparent trade off between manifold structure learning and clusters.**\n\nAs discussed in our general response, we want to emphasize that the distance preservation is not at odds with the cluster structure preservation. For clusters defined by the geodesic on a manifold, methods like HeatGeo, which are aimed at leveraging the underlying geometry of the data, should do better. In particular, as the reviewer points out, our method would remain faithful to the relationships between different clusters, in contrast to methods like t-SNE or UMAP, that may distort or inflate the separation. This is further illustrated on EMD interpolation experiment (Table 3) that aims at evaluating the quality of the relationship between different regions of the data. Experiments like the Gaussians on the swiss roll is a prominent example of clustering on a continuous manifold. In this case, the distance preservation and the cluster structure preservation are concordant objectives. \n\nThat said, some clustering tasks may indeed not perfectly fit the underlying manifold assumption, such as MNIST or COIL. We appreciate the reviewer’s suggestion of adding another high dimensional dataset, to put our method in perspective and we have added comparisons on MNIST and COIL20. You can see the full table in the PDF attached to the global response, this table will also be added to the main body of the text. On these datasets, t-SNE shows higher homogeneity and adjusted mutual information.\n\n**On the reviewer’s questions**\n\n>What is $V(x,sqrt(t))$ ?\n\nIt is the volume of a ball of radius $\\sqrt{t}$ centered at $x$, we added the definition.\n\n> Probably Vradhan could gain to be detailed more …\n\nIn the appendix C.3, we give a few examples on the application of the two-sided inequality and Varadhan’s formula, this would probably help the reader to understand. However we did not reference that section in the main body of the text, we have now fixed this. \n\n> the subsequent paragraphs in Preliminaries could be shortened or dropped…\n\nWe appreciate the reviewer’s advice on restructuring the preliminaries section and the related work. We believe we can shorten the subsection “Fast computation of Heat diffusion” and move some details to the appendix, and use the additional space to improve the fluidity of the related work section.\n\n> In section 4, if distance order preserving matters, why not compare to non-metric MDS too?\n\nWe added the comparison to non-metric MDS. The results are included in the PDF attached to the global response. Overall our method does better.\n\n> Can you clarify Section 4?...\n\nWe are working on improving the clarity of section 4. In particular, we added more intuition behind definition 4.1, and how it relates to the two-sided heat kernel inequality. We also fixed the notation on the norm $|\\cdot|$ vs $\\||\\cdot\\||_2$.\n\n>Section 5 is also wanting to embrace everything within limited space…\n\nSection 5 is a bit dense because it had to fit within the page limit. However in the appendix we provide the full details of the proof (Sec. A.2) and the unifying framework (Sec. A.1). In the appendix, before each proposition, we are adding a small paragraph explaining the results.\n\n>Ref 33: Cédric Villani has a twin brother?\n\nThank you for pointing this out, we fixed the reference. \n\nWe would like to thank the reviewer for their review of our paper. We believe we have answered all the great points raised by the reviewer in our author response. We respectfully ask the reviewer to reconsider their impression of the paper and potentially improve the given score if the raised concerns have been allayed. We thank the reviewer again for their time and we are also happy to answer any further questions that arise!\n"}, {"review_id": "6wox9RHjbP", "reviewer": "Reviewer_rWBe", "summary": "The authors propose heat geodesic embeddings, a novel approach to perform graph embeddings that is motivated by Riemannian geometry (Varadhan's formula and the parabolic Harnack inequality), which allows for directly linking heat diffusion on graphs to manifold distances. They performed experiments to demonstrate the efficacy of their approach. ", "questions": "- see sections above. ", "limitations": "- no negative societal impact. ", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- Originality: The proposed heat geodesic embedding and the approach is , to the best of my knowledge, original and novel. \n\n- Significance: It is certainly important to study spectral embeddings nowadays, and I find the authors' methodological contributions to be significant. \n\n- Quality: The quality of the theoretical derivations and experimental evaluations is sound, to the best of my knowledge. \n\n- Clarity: Despite the advanced mathematical material involved, the authors presented and motivated their approach in a clear and convincing manner. ", "weaknesses": "- there are no obvious weaknesses in the paper. The method is presented well and is straightforward. The authors performed experiments which covered most of the bases. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are pleased to read such an encouraging perspective on our work and thank the reviewer for their feedback. We remain at their disposal for any further clarification."}, {"review_id": "FgpqdQNMBh", "reviewer": "Reviewer_Yhnd", "summary": "In the paper \"A heat diffusion perspective on geodesic preserving dimensionality reduction\", the authors develop a new embedding method (HeatGeo), conceptually close to PHATE (Moon et al 2019) in that it performs MDS on some kind of diffusion distances. The authors make some strong claims about HeatGeo: (1) it is SOTA in preserving manifold distances; (2) it is SOTA in preserving cluster structures; (3) it can approximate both PHATE and t-SNE depending on parameters.\n\nI did not fully understand/follow all the math, but I found the paper nevertheless interesting as it gives a systematic/principled treatment of heat diffusion MDS methods. That said, I was not entirely convinced by the evaluations and especially claims (2) and (3) seem to me too strong and IMHO not sufficiently supported by the shown evidence. Also, many places were not really clear to me. Overall I am somewhat on the fence about this paper: it seems it is not fully ready yet. I am giving it a borderline score.", "questions": "See above", "limitations": "See above", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "* Systematic/principled mathematical treatment (which admittedly I did not fully grasp)\n* Benchmarks of several embedding methods on artificial datasets\n* Applications to real RNAseq data\n* Relationships between HeatGeo and existing methods (PHATE/t-SNE)", "weaknesses": "* Relationship to PHATE is not very clear\n* Very bold claims that are not fully supported by the evidence\n* Insufficient evaluation, in particular many embeddings are not shown\n* Runtime is not discussed\n\nNow in more detail:\n\n* In spirit, the method seems very similar to PHATE. First, kNN graph is used to construct NxN distance matrix between all points (using heat diffusion). Second, MDS is used to produce the embedding. Of course there are many small differences from PHATE in how exactly all of this is setup up, and the way this paper is written, it comes across as more mathematically principled than PHATE. But in the end, I was still confused about what is the MAIN conceptual/algorithmic difference. That's very important in order to judge the relative novelty of HeatGeo. Is it conceptually a minor tweak to / simplification of PHATE, or is it conceptually a very different thing?\n\n* One bold claim is that HeatGeo outperforms other methods in preserving clusters (supported by Table 2). This is conceptually surprising to me, because I expect some trade-off between preserving geodesic distances and preserving clusters (similar trade-off was discussed in Ref [2]: https://jmlr.org/papers/v23/21-0055.html), so it's odd if the same method (HeatGeo) wins in both competitions. Also, we know that methods like t-SNE and UMAP really excel in representing clusters in high-dim data (just try them on MNIST) but more geodesic/MDS methods like Isomap or PHATE perform worse (again, just try them on MNIST). \n\n  To be more convinced here, I would like to see the embeddings corresponding to Table 2. Also, I don't think \"Tree\" is a good dataset for clustering because we don't expect k-means to pick up the tree branches: just look at HeatGeo if Tree in Figure 1 -- How can k-means separate the branches? I don't really understand why HeatGeo performs well here in Table 2. Rather than the Tree, I would like to have some toy dataset with 100D Gaussians differently far away from each other. I would also love to see MNIST (a subset of MNIST).\n\n* Another bold claim is that HeatGeo approximates SNE/t-SNE. I was very surprised to read that, because I have never seen SNE/t-SNE formulated as an MDS problem, and assumed it is impossible. However, I don't really see much evidence that this is true. Can you do HeatGeo of MNIST that would look like t-SNE? In the EB dataset (Figure 3), can you change HeatGeo params to make it look like t-SNE? The same for the PBMC dataset (Figure 4)? \n\n* Figure 1: This is the main motivational figure showing that HeatGeo is very good at unrolling Swiss roll, whereas t-SNE creates \"artifical clusters\". Personally, I think it does not make much sense to apply t-SNE to an intrinsically 2D dataset because t-SNE is designed for reducing dimensionality (that's why it uses t-distribution kernel, unlike SNE). When working with intrinsically 2D data, it makes much more sense to use SNE instead of t-SNE. And indeed, SNE is much better at unrolling Swiss roll than t-SNE (example: https://jlmelville.github.io/smallvis/swisssne.html). As the authors discuss SNE in their paper anyway, I think it would be good to add SNE to the comparison tables and also to Figure 1 (can be done e.g. via openTSNE). By the way, Isomap can also unroll the Swiss roll perfectly.\n\n* As the HeatGeo relies on MDS, it constructs a dense NxN distance matrix and so has O(N^2) complexity. All the datasets used in the paper are very small. This is unlike t-SNE/UMAP which can be run on tens of millions of points easily. I think this has to be discussed, e.g. in Section 7.\n\n\nSMALLER POINTS AND QUESTIONS\n\n* line 124: classic MDS does not optimize objective (5).\n\n* line 148: Ref [2] shows that UMAP's effective loss function is very different from its stated loss function, so this formula is actually misleading. See also https://proceedings.neurips.cc/paper/2021/hash/2de5d16682c3c35007e4e92982f1a2ba-Abstract.html and https://openreview.net/forum?id=B8a1FcY0vi which develop this argument.\n\n* line 194: \"small t can lead to disconnected graphs\" -- I don't understand this sentence, I thought the graph does not depend on t, no? You first construct the graph and the laplacian L, and only later you use the diffusion time t.\n\n* line 211: Algorithm 1 -- I don't really understand how L is constructed. What is the underlying graph? Is it a kNN graph, what are the weights? What is the k value (number of neighbors)? I think this is not really spelled out in the text. PHATE used some complicated setup with alpha-kernels, are you also using this here?\n\n* Table 1: this table is missing t-SNE and UMAP rows. Also, I would like to see SNE in there, and Isomap.\n\n* Figure 3 caption: \"UMAP and t-SNE do not capture the continuous manifold\" -- actually UMAP looks pretty similar to HeatGeo to me. Also, Ref [2] argues that higher attraction in t-SNE leads to better manifold preservation, so UMAP is better at showing the manifold (which is what we see in your Fig 3) but one can do even better by using higher exaggeration in t-SNE. Maybe worth mentioning.\n\n* line 314: \"HeatGeo interpolated to PHATE for \\rho->1\" -- I was confused by that. What is rho? It wasn't really discussed before, I don't see it in section 5 that discusses relation to PHATE. I do see \\rho in Algorithm 1 but not in the text of Section 4. So it's mysterious.\n\n* Figure 4 does not show t-SNE, but Figure 5 in the Appendix does. I am surprised by how bad t-SNE looks. I have seen t-SNE of various PBMC RNAseq datasets many times, and t-SNE is usually very good at separating cell types. But in your appendix figure t-SNE looks very poor. Why is that?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their thorough analysis of our manuscript. In the following we will address the main points raised by the reviewers. Typos or undefined notation were fixed.\n\n**Relationship to PHATE is not very clear.** Our method follows the same main steps as Isomap and PHATE; 1) Graph construction, 2) Distance computation 3) MDS. Hence algorithmically, these methods are similar. Amongst the differences with PHATE is the use of the continuous dynamic (the heat kernel), the volume regularization, decoupling of the denoising. As you mentioned a key contribution is connecting diffusion-based distances to heat flow on the manifold, providing explicit bounds on the manifold geodesic distance. Thus we showcase a template that uses the framework of ISOMAP and PHATE but utilizes heat diffusion directly. \n\n**HeatGeo and clustering.** We would like to point the reviewer to the global response where we address the HeatGeo and clustering. (1) We believe that distance preservation is not at odds with the cluster structure preservation; (2) we used a mixture of Gaussian on the swiss roll (a combination of manifold and cluster); (3) we added evaluation on COIL-20 and MNIST. \n\n>Figure 1: This is the main motivational figure showing that HeatGeo is very good at unrolling Swiss roll…\n\n**On Figure 1.**  Thank you for providing us with this link, indeed it is a great example, however the results using 3k observations are not really convincing. In Fig. 1, we wanted to present results from the methods that are the most used in practice (although diffusion maps is not the most popular, we added it since it is another diffusion-based method). Further, for Table 2, we used different hyper-parameters specific to each method and selected them according to their performance on a validation set. We also used multiple repetitions to assess the variability of performance, while Fig.1 only shows one experiment. That said, Fig. 1 also clearly shows that UMAP and t-SNE don’t perform well on the tree dataset. The orange class is split on UMAP and the blue class is heavily distorted on the t-SNE representation. To compare the methods, we think that the quantitative results are much more powerful, in particular Tab.1 on distance preservation of the manifold. We have updated our comparisons, and now provide quantitative results for Isomap and metric-MDS. More details about the experiment setup is available at the beginning of Section 6, and the code is available [here](https://anonymous.4open.science/r/anon-heatgeo-CE2A/cluster_pred.py). \n\n**HeatGeo approximates SNE/t-SNE.** We would like to emphasize that our method can be reconfigured to be similar to a modification of SNE. We refer to proposition 4.3 for the details on the configuration. It is not true that HeatGeo directly approximates SNE, we will review the manuscript to make sure this statement is clear. We think the similarity between the two methods is clearer from a theoretical perspective. If we have space, we will move parts of the proof of prop. 4.3 to the main body of the text. The main point of this proof is that because the heat-geodesic dissimilarly uses the logarithm we can frame it as minimizing a KL loss (similar to SNE). The proof is provided in the supplementary material in section A.2.\n\n**On the graph construction.** The Laplacian is constructed from an affinity kernel. It is clearer in the full Alg.2 provided in the supplement, we will clarify it in the main body of the text. In our experiments we used the alpha-decay kernel to construct the graph. We experimented with other kernels and did not see a significant difference in accuracy. For example, in Fig. 13 we report the accuracy on the swiss roll for 4 types of kernels. \n\n**Complexity of HeatGeo.** Please refer to the global response for our comment on the algorithm’s complexity.\n\n> line 124: classic MDS does not optimize objective (5).\n\nYou are correct, we use the metric MDS problem defined in [1], we made it clearer in the manuscript. \n\n>Table 1: this table is missing t-SNE and UMAP rows. Also, I would like to see SNE in there, and Isomap.\n\nIn Tab.1 we compare the distance matrix used by a method against the ground truth distance matrix. UMAP, SNE, and t-SNE, do not rely on a distance matrix, but rather on an affinity kernel. So we cannot compare with them. We did a comparison against the shortest path distance used in Isomap, we noted it “Shortest path”. In this table, the focus is on the distances, not the embedding. \n\n>line 148: Ref [2] shows that UMAP's effective loss function…\n\nWe modified our statement and added the provided citations. \n\n>line 194: \"small t can lead to disconnected graphs\"...\n\nThis is a typo. The number of neighbors could cause a disconnected graph. We fixed it.\n\n> Figure 3 caption: \"UMAP and t-SNE do not capture the continuous manifold\"...\n\nWe agree that there is not a visually sharp distinction between the two. We modified the caption of Fig.3 accordingly. High exaggeration in t-SNE is indeed interesting! We have added a footnote to this effect. \n\n>I am surprised by how bad t-SNE looks...\n\nFor the clustering experiment, we try perplexity parameters of ​​10,30,100, and the figure is with perplexity 40, we are happy to try other parameters. But we did not observe that tSNE is good at separating this PBMC dataset. We believe there exist parameters where t-SNE performs well on PBMC, but it may take more work to find those. \n\nWe thank the reviewer for their valuable feedback and great questions. We hope that our rebuttal fully addresses all their important salient points, and kindly ask the reviewer to potentially upgrade their score if the reviewer is satisfied with our responses. We are also happy to answer any further questions that arise.\n\n[1] Zheng, J. X. et al. (2018). Graph drawing by stochastic gradient descent. TVCG\n\n[2] Moon, K. R. et al. (2019). Visualizing structure and transitions in high-dimensional biological data. Nature biotechnology."}]

Author rebuttal:
We would like to thank the reviewers for their time and valuable feedback when reviewing our paper. These suggestions led us to make a few modifications, which we believe greatly improved the quality of the manuscript. We now respond below to the main clarification points grouped by theme.

- **On our key contributions in unification of embeddings based on the idea of heat diffusion** We wish to highlight that one of our key contributions is to examine and interpret manifold-preserving embeddings based on properties of heat diffusion on the manifold. The key results we use to do this from geometry are Varadhan’s formula and Harnack’s inequality (which gives a two-sided bound on manifold distances based on heat diffusion). In particular, this perspective provides a better understanding of PHATE, diffusion maps, and tSNE. Then, we propose an embedding based on the heat kernel that makes this link more direct, i.e. the methods PHATE, diffusion maps and tSNE all use a random walk operator that indirectly models heat diffusion in discrete time whereas the heat kernel is a more direct continuous embodiment. As noted by reviewer Yhnd, this indeed makes our embedding “somewhat similar” (on purpose) to other manifold-preserving embeddings as it is a unification of the theory in those methods. 


- **On distance preservation and clusters**
The reviewers Yhnd and QP8B were skeptical about the clustering evaluation. In this manuscript we are proposing and unifying manifold-distance preserving embeddings. We want to note that distance preservation is in no way at odds with cluster structure preservation. Indeed, clusters are formed by data points that are close (low manifold distances) and separations between clusters are indicated by high distance. In fact, the distance preserving methods also keep the relationships as our method correlates well with the intrinsic distance on the manifold (see Tab.1). Indeed, a method that would perfectly recover intrinsic distances should also recover the resulting clusters perfectly. As such, for clusters defined by the geodesic distance on a manifold, methods like HeatGeo should do better, and indeed PHATE [2] showed superior clustering performance as compared to other methods like UMAP/tSNE which do not preserve distances. This is further illustrated in our EMD interpolation experiment (Table 3) that aims at evaluating the quality of the relationship between different regions of the data, and in Table 1 that reports the correlation between the different distances and the geodesic.
In our clustering experiments, we test a single cell PBMC dataset (known to have clusters), a mixture of Gaussian distributions on the swiss roll, and a tree with five branches. We have added experiments on MNIST and COIL-20; two datasets for which the manifold hypothesis may not be respected. Our methods show greater accuracy on the swiss roll and tree datasets, and competitive performance on the other datasets (see the full results in Tab.1 of the attached PDF). 

- **Additional comparisons to ISOMAP and non-metric MDS.** The reviewers Yhnd and QP8B mentioned the lack of comparison with Isomap, and reviewer QP8B also suggested using non-metric MDS. We followed their advice and added quantitative comparisons with non-metric MDS, and Isomap. For the MDS step in Isomap, we used the same optimization algorithm as our method, in order to accurately compare the two algorithms. Overall our method HeatGeo performs better than non-metric MDS and Isomap (see the preliminary results in the attached PDF). We will add visualizations for these methods in the appendix. This is to be expected in noisy data since shortest path distances (which Isomap uses) are highly sensitive to noise, and non-metric MDS ordinally preserves Euclidean and not manifold distances. 
	
- **On computational complexity.** 
Reviewer Yhnd asked for a discussion on complexity. We do mention the complexity of the approximation of the heat kernel (line 114), and we will add a discussion on the total complexity in the manuscript. While the most costly step of our current implementation constructs a full distance matrix for quadratic $O(n^2)$ cost, this is not strictly necessary. Similar to other methods, if we are willing to trade off accuracy for speed, we can truncate at some radius around each point to get a sparse distance matrix. If we do so then our algorithm will be $O(E)$ where E is the number of distances smaller than the chosen radius. This would give the same complexity (both time and space) as UMAP and TSNE algorithms, a log-linear-time sparse graph construction, a log-linear Chebyshev diffusion, followed by a (difficult to analyze) linear SGD optimization. We note that we use the stochastic gradient descent metric-MDS [1], which is much faster than using e.g. the SMACOF algorithm. In general our methods can also be made faster using landmarking [2, 3]. We also ran additional experiments to compare the computation time of the different methods. The preliminary results can be found in Tab. 3 of the attached pdf, which shows that our method is the fastest when using Chebyshev polynomials. We will show additional results by varying the number of datapoints.

- **We improved the readability and structure of the manuscript.** We did a thorough revision on the manuscript and corrected typos or misprints, improving the overall clarity of the manuscript.

[1] Zheng, J. X., Pawar, S., & Goodman, D. F. (2018). Graph drawing by stochastic gradient descent. IEEE transactions on visualization and computer graphics, 25(9), 2738-2748.

[2] Moon, K. R., van Dijk, D., Wang, Z., Gigante, S., Burkhardt, D. B., Chen, W. S., ... & Krishnaswamy, S. (2019). Visualizing structure and transitions in high-dimensional biological data. Nature biotechnology, 37(12), 1482-1492.

[3] Silva, V., & Tenenbaum, J. (2002). Global versus local methods in nonlinear dimensionality reduction. Advances in neural information processing systems, 15.

Meta-review: Reviews on this paper are mixed but tending positive; the most positive and most negative reviews both have some aspects that make them less than reliable. 

In the end, the AC decided to suggest accepting this paper.  Although the approach here is clearly related to many others in the diffusion maps literature, the idea certainly caught the interest of all the reviewers and generated discussion.  Moreover, the mathematical formulation here is different than motivation for past diffusion maps constructions of embeddings.
Decision: Accept (poster)