Learning Distances from Data with Normalizing Flows and Score Matching
Authors: Peter Sorrenson
Daniel Behrend-Uriarte
Christoph Schnoerr
Ullrich Koethe
Venue: ICML 2025 poster (2025)
Keywords: density-based distance
Fermat distance
Riemannian geometry
representation learning
normalizing flows
score matching
Primary area: general_machine_learning->representation_learning
TL;DR: Our work improves estimation of Fermat distances by combining normalizing flows, score-based models, geodesic smoothing, and a new dimension-adapted Fermat distance for better scalability.

Abstract: Density-based distances (DBDs) provide a principled approach to metric learning by defining distances in terms of the underlying data distribution. By employing a Riemannian metric that increases in regions of low probability density, shortest paths naturally follow the data manifold. Fermat distances, a specific type of DBD, have attractive properties, but existing estimators based on nearest neighbor graphs suffer from poor convergence due to inaccurate density estimates. Moreover, graph-based methods scale poorly to high dimensions, as the proposed geodesics are often insufficiently smooth. We address these challenges in two key ways. First, we learn densities using normalizing flows. Second, we refine geodesics through relaxation, guided by a learned score model. Additionally, we introduce a dimension-adapted Fermat distance that scales intuitively to high dimensions and improves numerical stability. Our work paves the way for the practical use of density-based distances, especially in high-dimensional spaces.

Learning Distances from Data with Normalizing Flows and Score Matching
Peter Sorrenson' Daniel Behrend-Uriare
Christoph Schnörr
Ullrich Köthe'
Abstract
100
density-weightee shortest path
Density-based distances (DBDs) provide a prin- cipled approach metric learning by defining distances in terms of the underlying data dis- tribution. By employing Riemannian metric that increases in regions of low probability den- sity, shortest paths naturally follow the data man- ifold. Fermat distances, a specific type of DBD, have attractive properties, but existing estimators based nearest neighbor graphs suffer from poor convergence due to inaccurate density estimates. Moreover, graph-based methods scale poorly to high dimensions, as the proposed geodesics are often insufficiently smooth. We address these challenges in two key ways. First, we learn den- sities using normalizing flows. Second, re- fine geodesies through relaxation, guided by a learned score model. Additionally, introduce dimension-adapter Fermat distance that scales intuitively to high dimensions and improves nu- merical stability. Our work paves the way for the practical use of density-baat based distances, espe- cially in high-dimensional spaces.
True gendesic
0.75
0.50
0.25
-0.25
-0.75
-1.00
-0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
1.25
Figure Our method Iearns distance metrics that respect un- derlying data distribution. Here, we show shortest paths between two points complex Gaussian mixture distribution. The ground truth geodesic (dashed red) and our graph-based approximation (solid blue) both naturally follow high-density regions where data samples (green dots) concentrated, rather than taking direct Euclidean path. This behavior is desirable for many machine learn- ing tasks, captures the intrinsic structure of the data manifold.
1. Introduction
Metric learning is fundamental to many machine learning tasks, from clustering to classification. The goal is to deter- mine a distance metric that accurately measures the similar- ity or dissimilarity between data points. While traditional approaches map data points into spaces where Euclidean dis- tance can be directly applied, these methods are inherently limited by the geometric constraints of Euclidean space, making impossible to represent arbitrary distance relation- ships between sets of points.
involves defining a Riemannian metric the data space and solving for geodesic distances. This allows for arbitrary distance relationships, but it requires choosing appropri- ate Riemannian metric. Fermat distances (Groisman et al., 2022), which are type of density-based distance (Bousquet et al., 2003), offer an elegant solution: the metric should conformal (a multiple the identity matrix) and in- versely proportional the probability density. Intuitively, this means that distances are stretched in regions where data is sparse and compressed where data is dense. In this way, geodesics pass through high-density regions of the data, respecting any inherent manifold structure (see Figure 1), This approach is consistent with Fermat's 'ss principle least time in optics, where the inverse density functions like a refractive index.
more flexible but computationally expensive approach
'Computer Vision and Learning Lab, Heidelberg University, Germany "Image and Pattern Analysis Group, Heidelberg Uni- versity, Germany. Correspondence to: Peter Sorrenson <pe- ter.sorrensonee @ @mmiil.ommm......
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Learning Distances from Data with Normalizing Flows and Score Matching
This crnnsty--bsed approach to metric learning has evolved significantly over the past two decades. The foundational work of Bousquet et al. (2003) introduced the concept of conformal density- based metrics, though their focus was on semi-supervised learning rather than computing geodesics. Sajama & Orlitsky (2005) took the first steps toward prac- tical implementation by using kernel density estimation and nearest-neighhot graphs for geodesic computation. A breakthrough came from Bijral et al. (2011) 1), who approxi- mated the density as inversely proportional to a power of the Fuclidean distance between close neighbors. This elegant approximation sparked numerous theoretical extensions and practical implementations (Chu et al., 2020; Little et al., 2022; Groisman et al., 2022; Moscovich et al., 2017: Alam- gir & Von Luxburg. 2012: Mckenzie & Damelin, 2019; Little et al., 2020), becoming known as power-weighted shortest path distances (PWSPD) (Little et al., 2022) Fer- mat distances (Groisman et al., 2022; Trillos et al., 2023). While Hwang et al. (2016) provided theoretical convergence guarantees for these methods, notably absent from all previ- ous work is any comparison to ground truth distances, even in simple cases where they can be computed exactly. Our empirical analysis reveals that existing approaches exhibit poor convergence in practice, highlighting a significant gap between theoretical guarantees and practical performance.
densities and their gradients (scores).
In summary, our contributions are fourfold: First, we in- troduce a numerically stable relaxation method to compute ground truth geodesics and distances. Second, we intro- duce the use of normalizing flows to obtain more accurate edge weights for estimating Fermat distances, achieving sig- nificantly faster convergence rates than previous methods. Third, we overcome the limitations of graph-based methods in high-dimensional spaces through novel smooth relax- ation scheme using score matching. Finally, we propose a Fermat distance that naturally scales with dimensional- ity, making density- distances practical for real-world, high- digh-dimensional applications. Our approach bridges the gap between theoretical elegance and practical applicability in metric learning.
2. Background
2.1. Riemannian Geometry
Riemannian geometry provides the mathematical foundation for understanding curved spaces and measuring distances within them. While Euclidean geometry suffices for flat spaces, many real-world datasets lie on curved manifolds where Euclidean distances can be misleading. Riemannian geometry offers powerful tools for analyzing such spaces through the metric tensor, which can encode different ge- ometric structures ranging from positive curvature (like a sphere) to negative curvature (like hyperbolic space), or more complex combinations of both.
Other notable works which use density-based metrics to inform the learning of data manifolds include metric flow matching (Kapusniak et al., 2024) and the Riemannian au- toencoder (Diepeveen et al., 2024), which uses a metric defined in terms of score functions.
More formally, a Riemannian manifold M is a smooth n- dimensional space that locally resembles RR",, In this work, we make the simplifying assumption that M = R". As result the tangent space Tp. which is linear approximation of the manifold at p. is simply R" for all p. The metric tensor 9p : TR x TR" -> R is a smoothly varying positive- definite bilinear form on the tangent space. It defines an inner product in the tangent space, allowing us measure lengths and angles.
To establish ground truth geodesics and distances, we de- velop a numerically stable relaxation method that directly solves the geodesic differential equations. By carefully han- dling the reparameterization of paths to maintain constant Euclidean velocity and employing appropriate numerical op- timization techniques, we can compute accurate geodesics even in challenging cases with large density variations. This enables us, for the first time, to quantitatively evaluate the accuracy of different approximation methods by comparing them to true geodesics in distributions with known densities, such as mixtures of Gaussians.
For tangent vector E Tp the norm is defined as ||||| = V9p(t, v), providing a notion of length in the tan- gent space. Just as we measure the length of curve in Euclidean space by integrating its speed, the length of a smooth curve r(t) on the manifold is given by I (t)( dt, where (t) is the derivative of (t) with respect to f. This in- tegral accumulates the infinitesimal lengths along the curve, accounting for how the metric varies across the manifold.
With this ground truth framework place, we develop two key improvements to existing methods. First, by employing more accurate edge weights obtained through normalizing flows, we achieve much faster convergence in graph-based approximations. Second, we address the limitations of these graph-based methods in higher dimensions, where paths be- come increasingly rough due to the curse of dimensionality, by using a relaxation scheme to smooth trajectories. Inter- estingly, we find that training models via score matching is more effective than using normalizing flows for this purpose, revealing trade-off between the accuracy of learned log
The distance between two points on the manifold is defined as the minimal length of curve connecting the points. The curves that minimize this length are known as geodesics, which are the generalization of straight lines in Euclidean space to curved spaces. These geodesics satisfy a second-
2

---

Learning Distances from Data with Normalizing Flows and Score Matching
order differential equation known as the geodesic equation (see Appendix A. In our context, these geodesics will 05-
Uaodaa Jd
Sadid (2 1/A
represent optimal paths through the data space that respect the underlying probability density structure.
44-
2.2. Fermat Distances
Having established the framework of Riemannian geometry, we now turn to density-basee baseed distances (DBDs), which pro- vide a way to define distances between points sampled from a distribution by utilizing the underlying probability density function. This approach can be particularly useful for tasks such as clustering and pathfinding. If designed well, a DBD can capture desirable properties such as:
Figure 2. Geodesics of the standard normal distribution, projected into 2D, taken between two orthogonal points at distance / from the origin. The paths are scaled by vD to enable comparison. Left: Unscaled metric (same temperature for all dimensions) leads to sharply curved geodesics in high dimensions. Right: Scaled metric (temperature equal dimension) leads to consistent geodesics in all dimensions (note that the lines are overlapping).
1. Shortest paths between points should pass through sup- ported regions of the data. If there is a region with no data between two points, the shortest path should deviate around it to pass through regions with more data.
The length of a smooth curve (t) is therefore given by
L(>) - / 1G>=
I1( di p(T(1))
(2)
This integral accumulates the Euclidean length element [|((t)],, weighted inversely by the density along the path. As a result, segments of the path that pass through low- density regions contribute more to the total length than those through high-density regions.
2. Data should cluster according to its modes. If two data points belong to different modes with little data connecting them, the distance between them should be high.
Fermat Distances Which Scale with Dimension While the Fermat distance framework is theoretically elegant, set- ting 3 = (as is common in previous work) leads to unintu- itive behavior in higher dimensions. Consider the standard normal distribution in D dimensions. Due to rotational sym- metry, the geodesic connecting #1 and 222 will lie in the plane spanned by these two points and the origin. Without loss of generality, suppose this plane aligns with the first two axis directions. The density restricted to this plane is the same as in the two-dimensional case, up to constant factor. Therefore, geodesics in the D-dimensional standard normal distribution are the same as in the two-dimensional standard normal distribution, once we account for the rotation into the plane of z1 and Z2.
One way to define a DBD that satisfies these properties is through Fermat distances. Named after Fermat's principle of least time in optics, these distances model the shortest paths between points by considering the inverse probability density (or a monotonic function of it) as analogous to the refractive index.
In the language of Riemannian geometry, we define a metric tensor of the form
(u, v) (u, = f(p(x))²*
(1)
where p the probability density, f is monotonic function, and (,,)) is the Euclidean inner product. For example, in mixture of two Gaussians, this metric would make paths that cross the low-density region between the modes longer than paths that go through high-densit regions, effectively capturing the natural clustering in the data. This formu- lation ensures that regions with higher probability density have lower "refractive index," " guiding shortest paths through denser regions of the data. In practice, is convenient to choose {p) = pº, where 3 is an inverse temperature param- eter that controls the influence of the density on the metric. Since our metric enacts a simple scaling of the Euclidean inner product, it is a conformal metric, which locally pre- serves angles and simplifies key geometric objects such as the geodesic equation (see Lemma A.1 and Theorem A.2).
However, despite this geometric similarity, there's a cru- cial scaling issue: the typical distance from the origin of points sampled from the D-dimensional standard normal distribution is much larger, with a mean of VD. As a result, geodesics will start and end in exponentially low-density regions relative to the two-dimensional standard normal. This forces most of the trajectory to make a dramatic detour towards and then away from the high-density region near the origin, creating increasingly extreme curved paths as dimension grows.
To address this issue, if we scale the temperature linearly with dimension (3 = 1/D), the results become consis- tent with the wo-dimensional case (see Figure 2). In Ap-
3

---

Learning Distances from Data with Normalizing Flows and Score Matching
Semee Derstr
1-99 Crach
Figure 4. Illustration of the power weighted graph approximation. Left: Sample points (dots) and density (gradient) from 2D stan- dard normal. Center: k-nearest neighbor connections between points, with edge weights determined by cubic Euclidean dis- tances. Right: A shortest path (solid bluc) found by Dijkstra's algorithm compared the true geodesic (dashed red). Note how the discrete path approximates the continuous geodesic by follow- ing high-densit regions, but is not very smooth.
Figure 3. Comparison of geodesic parameterizations. Left: Con- stant metric speed leads to varying Euclidean step sizes, making numerical relaxation challenging. Right: Constant Euclidean speed enables uniform discretization and more stable numerics. Points along the paths indicate discretization steps. Note the slightly dif- ferent solutions: the left-hand path is approximately 0.1 1% longer.
pendix A.3, we explore this argument in more detail and generalize it to non non-Gaussian distributions. This scaling has two key benefits: First, it ensures that geodesic behavior remains intuitive across dimensions, with similar qualitative properties regardless of the ambient dimension. Second, it significantly improves numerical stability by reducing the effect of the probability density, leading more Euclidean- like geodesics. This effect is clearly visible in Figure 2, where the geodesics are much closer to straight lines in the scaled case.
To address these numerical challenges, we reparameterize the equation to maintain constant Euclidean speed instead of constant metric speed (see Theorem A.3). This simplifi- cation is particularly valuable for numerical methods like relaxation, where uniform spacing of points along the tra- jectory leads to better stability. Figure 3 illustrates how this reparameterization enables uniform discretization compared to the varying step sizes required by constant metric speed. The reparameterized equation differs only slightly from the original (losing just the factor of in the second term):
2.3. Geodesic Equation
F - 3(>(ip) G) + (=(()) = =0
(4)
From the metric, we can compute the geodesic equation, which characterizes the paths of minimum length between points. Just as a straight line in Euclidean space can be described by a second-order differential equation (i = 0), geodesics in our curved space satisfy a more complex equa- tion (see Theorem A.2). For geodesic trajectory (t), we have:
We stress that, while Equation (3) is easily obtained from standard results on conformal metrics (c.g.. Carroll (2004, Appendix G)), the reparameterized form in Equation (4) is, to our knowledge, novel contribution.
2.4. Weighted Graph Algorithm Based on Euclidean Distance
7- - 2/(((1)) : 4)4 + (6(1))11112 =
(3)
where s(x) plogp(z) is the score of the distribution (the gradient of the log density) and ||-I| is the ordinary Euclidean norm.
A practical approach computing Fermat distances is to discretize the problem using sample points. If we have large number of samples from p(z), we can approximate the density between close neighbors x1 and 22 constant and hence the shortest path between them as straight line. Following Bijral et al. (2011), the density is roughly pro- portional to the inverse of the Euclidean distance, raised to the power d, where d is the intrinsic dimension of the distribution:
&r
The geodesic equation has the property that the speed is constant with respect to the metric, meaning V 7) is constant. While mathematically elegant, this creates nu- merical challenges: when 7 passes from a high-probability region a cow-probability one, the Euclidean velocities can differ by orders of magnitude. For example, the probabil- ity density drops by a factor of 1000, the Euclidean velocity must increase by the same factor to maintain constant metric speed. This dramatic speed variation makes standard numer- ical integrators unstable and requires extremely small step sizes in density regions while potentially overshooting in high-density regions.
1
p(7)
(5)
1101
x22 dd
Hence the distance between Ti and 22 according to the metric is approximately proportional to power of the Eu- clidean distance
dist .22 x ||:1 x221 1d21
(6)

---

Learning Distances from Data with Normalizing Flows and Score Matching
This elegant relationship suggests simple algorithmic ap- proach: Bijral et al. (2011) proposes to approximate short- est paths between more distant points by constructing k-nearest neighbors graph of the samples with edge weights given by || - 22 || and applying Dijkstra's ss algorithm. Figure illustrates the method.
logllike training set:
f E= Pana(x)) po
(8)
2.6. Score Matching
While normalizing flows provide a way to estimate the den- sity directly, our geodesic equations actually depend on the score function s(x) V. log p(z). which is the gradient of the log- likelihoodd with respect to spatial inputs. Although one could obtain scores by differentiating the density model, we find that these derived scores are too noisy to be used effectively in practice. We speculate this reflects an inher- ent trade-off between optimal learning of the log density and its gradient, which we illustrate with simple examples in Appendix C.3. Therefore, we turn to score matching (Hyvärinen, 2005), which offers direct way estimate the score function from data, bypassing the need to estimate the density itself.
Hwang et al. (2016) gives consistency guarantees for this estimate, proving that a scaled version of converges to the ground truth distance with rate of exp(- (un (xp(-8)nn/(442)) ((34+2)), ) with n the sample size, d the intrinsic dimension, and Og a positive constant. Groisman et al. (2022) additionally show that the shortest paths converge to geodesics in a large sample size limit.
While this approach is appealing in its simplicity, can be easily implemented using standard libraries, and comes with theoretical guarantees, the density estimate is too poor to provide accurate estimates of Fermat distances in practice. In the experimental section 4.1.1 we show that paths com- puted with this approach converge extremely slowly to the ground truth geodesics as sample size increases, even in very simple 2-dimensional distributions such as the standard normal.
This family of methods includes several variants such as sliced score matching (Song et al., 2020), denoising score matching (Vincent, 201 and diffusion models (Song et al., 2021). For our purposes, we use sliced score matching due to its simplicity and ability to estimate the true score, rather than noisy score in denoising score matching diffusion methods.
This significant gap between theoretical guarantees and prac- tical performance leads us to speculate that the unknown constant 0g is very small, resulting in extremely slow con- vergence despite the exponential form of the bound. This observation motivates our development of more accurate density estimation methods in the following sections. We leave detailed theoretical investigation of this convergence behavior to future work.
Sliced Score Matching We can learn the score of p(x) by minimizing the score matching objective:
(9)
In high dimensions the trace term is expensive evaluate exactly, so sliced score matching uses a stochastic approxi- mation (the Hutchinson trace estimator (Hutchinson, 1989)):
2.5. Normalizing Flows
Normalizing flows (Rezende & Mohamed, 2015; Kobyzev et al., 2020; Papamakarios al., 2021) provide a powerful framework for learning complex probability distributions. The key idea is to transform simple base distribution (like standard normal) into the target distribution through series of invertible transformations. Formally, this is achieved via a parameterized diffeomorphism fo which maps from data to latent space. The probability density is obtained through the change of variables formula:
Bp(r)r(o) [e" V.so(z)) + $ (1=)
(10)
where is sampled from an appropriate distribution, typi- cally standard normal.
3. Methods
3.1. Ground Truth Distances through Relaxation
0.f ax
p=(x)= = T (fo(z)) det
(7)
We can solve Equation (4) in several ways to find geodesics between two points. The shooting method starts at one point with an initial velocity guess and integrates the equation forward, iteratively adjusting the initial velocity to reach the target point. Alternatively, the relaxation method dis- cretizes the interval [0, 1] into equal time steps of size h = 1/n, writing po: = ((ihh) for the curve values at these points. Starting from an initial path (either a straight line between endpoints or a more informed guess), it gradually updates the intermediate points while keeping the endpoints
where π a simple latent distribution. A common strat- egy to design normalizing flows is through coupling blocks (Dinh et al., 2014; 2017; Durkan et al., 2019), where fo is the composition of a series of blocks, each with triangular Jacobian. This makes the determinant easy to calculate, and hence the above formula is tractable.
Normalizing flows are trained by minimizing the negative
5

---

Learning Distances from Data with Normalizing Flows and Score Matching
Path Comparison (Graph Baseh-
Path Comparison (Rdaxatioo)
Algorithm 1 One step of relaxation (iterate until conver-
Porer wedy sho test path
Scure model Tuur goodesic
Density mode!
gence)
Ground tutn couuity
The pesceaii
Input: trajectory y = (p:) score function s, inverse temperature /3
1 {I,.... n 1}
shuffle(I) for i in mnn/on do
= (P - pi-1)
w; = 2 {
(s(()) m:)os]
-1
9, (fit) + +P;-11 + wi
end for Return: 4
Figure 5. Illustration of geodesic approximations by graph-based methods (left) and relaxation methods (right) 200 points sam- pled from the 2D standard normal. Left:: the weighte shortest path method (Section 2.4) produces the roughest path. The shortest path from an inverse-density weighted graph using model density (Section 3.2) is more direct, but not quite close to the true geodesic as the equivalent method using the ground truth density. The model defined as Gaussian with the same mean and covariance the data. Right: relaxation (Algorithm I) on the model score function results a very close approximation to the true geodesie.
Algorithm 2 Construction of brnnstty-weigtedd graph
Input: data matrix X, density function p. inverse temper-
ature 3, neighbors k. segments S
G = KNN k) {Construct k-nearest neighbor graph)
for edge (1, m) in 9 do
Ar = (Xm m X1)/S
for i = to S do
1/21/2 = X1 + (i - 0.5)
fixed to better satisfy the differential equation. Each method has its advantages: shooting more efficient for exploring geodesics from starting point, while relaxation is more ro- bust for finding paths between specific endpoints, especially in complex density landscapes where shooting might miss the target. We opt to use relaxation.
end for
9im - Ar Σ 1/p(U/-1/2)
end for Return: G
Equation (4) leads to the relaxation scheme given in Algo- rithm 1, which updates each point along the curve based on its neighbors and the local score function. The ran- dom ordering of updates helps avoid oscillatory solutions and improves convergence. See Appendix A for detailed derivation.
We use this distance measure to construct a k-nearest neigh- bors graph, with edge weights given by the approximate geodesic distances rather than Euclidean distances (see Al- gorithm 2). The resulting sparse graph structure allows us to efficiently compute distances between any pair of points using Dijkstra's 's shortest path algorithm. Due to the large dynamic range of the density values involved, we perform all distance calculations and graph operations in log space to maintain numerical stability.
3.2.
Graph Algorithm
To improve upon previous graph-based methods, replace their density estimates with those from normalizing flow trained on the data. For nearby points, we can approximate the geodesic connecting them as a straight line and estimate its length through numerical integration. Specifically, given two points Ii and x2,, we estimate their distance as:
3.3. Relaxation with a Learned Score Model
The relaxation algorithm described in Section 3.1 requires only the score function s(x), not the density itself. This observation allows us to use a score model trained directly via score matching, bypassing potential issues with density (11) estimation. In practice, we find this approach more effective than computing scores by differentiating normalizing flow (see experimental results in Section 4).
dist(r1, Re
1=1 po (1/1-1/22)
| 122l
1
(12)
S
pe 1/
To ensure robust convergence, we initialize the path using the graph-based method from the previous section. This initialization provides a reasonable approximation of the geodesic. helping the relaxation process avoid local minima and converge more quickly to the optimal path.
where { are equally spaced points along the line from I1 to 22 (i.e., yo = ω1 and ys = 22). and 1/-1/2 denotes the midpoint between consecutive points.
6

---

Learning Distances from Data with Normalizing Flows and Score Matching
Men Log Puth
70m
1071
Danwtyy Normeliinng Flow
Bldkklin Plwee
10"4
10"4
30
10F
10P
11P
32
10"
10P
10P
Saupls
Sauppk
Chrcle
Spinl
Spicali
Mas Log Fath Ratio
1
10-1
10"1
1074
10r4
10-4
Sumple
Semgle
Figure 6. The length of paths in learned and ground-truth density weighted graphs (orange and blue) converge to the ground truth at roughly the same rate with increasing sample size, while power weighttd graph method (green) converges a far inferior rate.
4. Experiments
4. 4.1.1. POOR CONVERGENCE OF POWER-WEIGHTSS METHODS
We have published the code to reproduce these experiments in the following GitHub repository: https: / //g/gttttttt github. com/viislea.m//eeratttt stannce.
We first examine the power-weighted graph method of Bijral et (2011), which constructs paths on k-nearest neighbor graphs. Despite theoretical guarantees, we find that perfor- mance improves extremely slowly with increasing sample size across all datasets (Figure 6, green lines).
Performance Metric To evaluate the quality of computed paths, we introduce the Log Path Ratio (LPR) metric that can be applied to any path & connecting two points I and We compute the path length ((p) using the ground truth density and compare it to the true geodesic distance:
4.1.2. DENSITY ESTIMATION IS THE BOTTLENECK
To identify the source of this poor convergence, we de- veloped several variant methods that use different nearest- neighbor density estimators while maintaining the same graph structure (see supplementary material for details). All variants show similarly poor performance, while their coun- terparts using ground truth densities converge rapidly. This clearly identifies density estimation as the primary bottle- neck.
L(p)
LPR(()) = log
(13)
dist(r1, 22)
This metric is zero only for the true geodesic and positive otherwise, providing principled way evaluate different methods when ground truth densities are available.
4.1.3. IMPROVED CONVERGENCE THROUGH BETTER DENSITY ESTIMATION
Having identified poor density estimation the is- sue, we replace nearest-neighbor estimates with learned normalizing flow model. Using this density estimate in Al- gorithm 2, we find that the resulting paths converge at nearly the same rate as those computed using ground truth densi- ties (Figure 6, orange and blue lines). This demonstrates that modern density estimation techniques can effectively close the gap between theoretical guarantees and practical performance. Full experimental details can be found in the
4.1. Convergence Analysis 2D Datasets
To evaluate convergence behavior, we analyze five two- dimensional datasets with known ground truth densities, ranging from simple unimodal distributions to complex mul- timodal mixtures (shown in Figure 6). For each dataset, we compute the average LPR over 1000 randomly sampled paths. Full experimental details and dataset descriptions are provided in the supplementary material.

---

Learning Distances from Data with Normalizing Flows and Score Matching
Unerled
Salee
mance, maintaining good convergence across dimensions for both scaled and unscaled metrics, though the scaled ver- sion (3 D) shows more reliable behavior (Figure 7), Detailed experimental results and analysis can be found in the supplementary material.
4.3. MNIST
We also perform preliminary experiments on MNIST within the latent space of an autoencoder, as described in Ap- pendix C.I. We find that we obtain interpretable and reason- able geodesics connecting digits within this space, and that digits cluster together in an unsupervised way.
Powcs
Nemndingg Ean
The Daudtz
Sonc Hold Rdanstcs
Figure 7. Mean LPR sample size 200.000. Graph-based based methods (first three) show rapidly degrading performance as di- mension increases, even with the scaled metric (3 = 1/D). In contrast, relaxation with a learned score function (red line) main- tains good performance across dimensions, effectively avoiding the curse of dimensionality.
5. Conclusion
We have presented a systematic study of geodesic compu- tation in density-based metrics, using distributions with known ground truth to rigorously evaluate different ap- proaches. Our investigation reveals several key insights:
supplementary material.
First, demonstrate that previous graph-based methods, despite their theoretical guarantees, converge impractically slowly even for simple distributions like the standard nor- mal. Second, we show that this poor convergence stems primarily from inadequate density estimation rather than the graph approximation itself. Third, we establish that modern deep learning approcche-nnomaalizinn flows for density estimation and score matching for gradients-can dramatically improve performance.
4.2. Scaling the Dimension
Machine learning datasets often lie in high-dimensional spaces, making dimensional scaling behavior crucial for practical applications. We investigate this by extending one of our two-dimensional datasets (the standard normal) up to 25 dimensions, testing both the scaled (3 = D/D) and unscaled (3 = I) versions of the metric.
4.2.1. GRAPII-BSEDD METHOΟS FAIL IN HIGHER DIMENSIONS
However, significant challenges remain. While our meth- ods work well on toy distributions with known densities, extending them to complex real-world data requires further research. Key directions for future work include:
All graph-based methods show rapidly degrading perfor- mance dimension increases, regardless of whether we use scaled or unscaled metrics (Figure 7). This failure is a direct consequence of the curse of dimensionality: with fixed sample size but increasing dimension, the probability of finding sample points near the true geodesic becomes vanishingly small. As result, graph-based paths become increasingly jagged approximations of the smooth ground truth trajectory.
Unifying normalizing flows and score models into single framework to improve efficiency
Developing theoretical understanding of why meth- ods converge effectively while previous approaches with consistency guarantees do not
Extending these methods higher-dimensiin densi- ties typical in real applications
4.2.2. SCORE- BASED RELAXATION OVERCOMES DIMENSIONAL SCALING
This work establishes a foundation for practical density- based distances, bringing theoretical elegance closer to prac- tical applicability in metric learning.
To address this fundamental limitation of graph-based meth- ods, we apply our relaxation algorithm (Algorithm 1) with learned score function. While one might expect to obtain scores by differentiating trained normalizing flow, find these derived scores too noisy for reliable convergence (see Appendix C.3 for analysis of this phenomenon). Instead, direct score estimation via sliced score matching provides smoother gradients that enable stable relaxation.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences our work, none which we feel must be specifically highlighted here.
This score-based approach dramatically improves perfor-
8

---

Learning Distances from Data with Normalizing Flows and Score Matching
Acknowledgments
Harris, C. R., Millman, K.J., Van Der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., TIylor, J., Berg, S., Smith, N. J., et al. Array programming with
This work is supported by Deutsche Forschungsgemein- schaft (DFG, German Research Foundation) under Ger- many's Excellence Strategy EXC-2181/1 390900948 (the Heidelberg STRUCTURES Cluster of Excellence). The authors acknowledge support by the state of Baden- Württemberg through bwHPC and the German Research Foundation (DFG) through grant INST 35/1597-1 FUGG,
NumPy. Nature, 585(7825)3357336 2020.
Hunter, J. D. Matplotlib: A 2d graphics environment. Com-
puting in science & engineering, 9(03)990-95, 2007.
Hutchinson, M. F. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Com- munications in MlttsticssSimullati and Computation, 186))1059-1076, 1989.
References
Alamgir, M. and Von Luxburg, U. Shortest path distance in random k-nearest neighbor graphs. In Proceedings of the 29th International Coference on International Conference on Machine Learning, pp. 1251-1258 2012.
Hwang, S. J., Damelin, S. B., and Hero III, A. O. Shortest path through random points. Annals of Applied Probabil- ity, 2016.
Hyvärinen, A. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
Ardizzone, I.., Bungert, T., Draxler, F., Köthe, U., Kruse,
J., Schmier, R., and Sorrenson, P. Framework for Eas- ily Invertible Architectures (FrEIA), 2018-2022. URL https: / github.. wttpm/iiiseeae/rr.......
Kapusniak, K., Potaptchik, P., Reu, T., Zhang, L.,, Tong, A., Bronstein, M., Bose, J., and Di Giovanni. F. Metric flow matching for smooth interpolations on the data manifold. Advances in Neural Information Processing Systems, 37: 13501 1-35042,, 2024.
Bijral, A. S., Ratliff, N., and Srebro, N. Semi-supervised learning with density based distances. In Proceedings of the Seenty-Seventh Conference on Uncertainty in Artifi- cial Intelligence, pp. 43-50, 2011.
Kobyzev, I., Prince. S. J., and Brubaker, M. A. Normalizing flows: introduction and review of current methods. IEEE transactions on pattern analysis and machine intel- ligence,
Bousquet, O., Chapelle, O., and Hein, M. Measure based regularization. Advances in Neural Information Process- ing Systems, 16, 2003.
2020.
Carroll, S. M. Spacetime and Geometry. Addison Wesley, 2004.
LeCun, Y., Bottou, L.,, Bengio, Y., and Haffner, Gradient- based learning applied to document recognition. Proceed- ings of the IEEE, 86(1
Chu, T. Miller, G. L., and Sheehy, D. R. Exact computation of a manifold metric, via Lipschitz embeddings and short- est paths on a graph. In Proceedings of the Fourteenth Annual ACM-SIAN Symposium on Discrete Algorithms, pp. 411-425 SIAM, 2020.
2002.
Little, A., Maggioni, M., and Murphy, J. M. Path-based spee- tral clustering: Guarantees, robustness outliers, and fast algorithms. Journal of machine learning research, 21 16):1-66,7 2020.
Diepeveen, W., Batzolis, G., Shumaylov, Z., and Schönlieb, C.-B. Score-based pullback Riemannian geometry. arXiv preprint arXiv:2410.0/95
Little, A., McKenzie, D.,, and Murphy, J. M. Balancing ge- ometry and density: Path distances high -iseeniinnsi data. SIAM Journal on Mathematics of Data Science, 4 (1):72-99, 2022.
2024.
Dinh, L., Krueger, D., and Bengio, Y. Nice: Non-liness independent components estimation. arXiv preprini arXiv: 2014.
Mckenzie, D. and Damelin, S. Power weighted shortest paths for clustering Euclidean data. Foundations of Data Science, 13)):307-327, 2019.
Dinh, I.,, Sohl-Dickstein, J., and Bengio, S. Density esti- mation using real nvp. In International Conference on Learning Representations, 2017.
Moscovich, A., Jaffe, A., and Boaz, N. Minimax-optimal semi-supervised regression on unknown manifolds. In Artificial Intelligence and Statistics, pp. 933-942. PMLR, 2017.
Durkan, C., Bekasov, A., Murray, I., and Papamakarios, G. Neural spline flows. Advances in neural information processing systems, 32, 2019.
Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., and Lakshminarayanan, B. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1 1457))1-64, 2021.
Groisman, P., Jonckheere, M., and Sapienza, F. Nonhomo- geneous Euclidean first- pest-passage percolation and distance learning. Bernoulli, 28(1):255-276 2022.

---

Learning Distances from Data with Normalizing Flows and Score Matching
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z.,, Gimelshein, N., Antiga, L., et al. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32, 2019.
Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al. Scikit-leam: Machine learning in Python. the Journal of machine Learning research, 12:2825-28333 28300, 2011.
Rezende, D. and Mohamed, S. Variational inference with normalizing flows. In International conference on ma- chine learming, pp. 1530-1538. PMLR, 2015.
Sajama and Orlitsky, A. Estimating and computing density based distance metrics. In Proceedings of the 22nd inter- national conference on Machine learning, pp- 760-767, 2005.
Song, Y.. Garg, S., Shi, J.. and Ermon, S. Sliced score matching: A scalable approach to density and score timation. In Uncertainty in Artificial Intelligence, pp. 574-584. PMLR, 2020.
Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Er- mon, S., and Poole, B. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations, 2021.
Trillos, N. G., Little, A., McKenzie, D., and Murphy, J. M. Fermat distances: Metric approximation, spectral con- vergence, and clustering algorithms. arXiv preprint arXiv:2307.05750, 2023.
Vincent, A connection between score matching and de- noising autoencoders. Neural compulation, 23(7):1661- 1674, 2011.
Virtanen, P., Gommers, R., Oliphant, T. E.,, Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P.,
Weckesser, W., Bright, J., et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature methods, 2613):261-272, 2020.
10

---

Learning Distances from Data with Normalizing Flows and Score Matching
A. Derivations
A.1. Geodesic Equation
We begin by deriving the geodesic equation for our conformal metric. This is well-known result in differential geometry (see e.g., Carroll (2004, Appendix G)). The key steps are computing the Christoffel symbols and then simplifying the resulting equation.
Lemma A.1 (Christoffel Symbols for Conformal Metrics), For a conformal metric of the form gej = X28.g, the Christoffel
symbols are:
alog X Or
θ Olog X
alog X Oz
03
5%
8
(14)
Oxi
Proof. The Christoffel symbols are given by:
Oggi Orch
0gi Oxl
Ogjk Ox!
jk = 29
+
(15)
For our metric:
1A2
Ox2
OX3
r
8t
8i:
Ox 8jk
(16)
212
Oz
Oxi
8t θ @log
Ø Olog X
Olog
dik
8sk
(17)
Ork
Ozi
Ox
The result follows by simplifying the Kronecker delta products.
Theorem A.2 (Geodesic Equation). For metric of the form Piz 1-p2--00pp deji where is a probability density and B3 is constant, the geodesic equation in vector notation is:
ल - -3((() 4) + BB(7)
(18)
where s V log p is the score function of the probability density.
Proof. For our metric, symbols:
- and hence
Blog Dn
-3
3.logg -Bs. Substituting into the formula for Christoffel 8x
2 BB - Bs,0g 0
(19)
The geodesic equation 7 +
jk 13gh then becomes:
74 - B Bs + B505 3 d
(20)
which simplifies to the desired vector form.
Theorem A.3 (Constant Euclidean Speed Parameterization). Under a reparameterization that maintains constant Euclidean speed, the geodesic equation becomes:
の 3- ())()))))+++++++ + Bs( =
(21)
Proof. Let f [0, 1| + [0, 1] be strictly increasing diffeomorphism and define vp(u) = (w)- The derivatives are
related by:
- 6f2 + 45 = =
(22)
(23)
For constant Euclidean speed, we require:
du (1
d
= -b- x 8=
(24)
Substituting into the geodesic equation and solving for j yields:
)= /2
(25)
The result follows by substituting back and simplifying.

---

Learning Distances from Data with Normalizing Flows and Score Matching
A.2. Relaxation Scheme
The central finite difference approximations for the first and second derivatives are
65 Pi+1
-Pi-1
(26)
2h
and
Pi+1
24p
+ 4ps-1
Rd
(27)
122
By substituting this into the differential equation (Equation (21)) we have
Pi+1
24, + pi-1 ((((pp))) (b.) + Bs (=(())))
(28)
h2
where b is estimated using finite differences. We can rearrange this to
1²8 Pi Pi+1 + Pi-1 + (s(r)) (s(ppt (())))))
(29)
and use the equation as an update rule, updating each position of the curve except the endpoints at each iteration.
Dividing by very small h could lead to numerical instability. We can avoid dividing and multiplying by by the following update rule. First define
P1+1 Pi-1
U:
(30)
2
then update with
Pit1 + Pi-1
B
(s(p))llull- (s(pp) v))s))
(31)
2 2
A.3. Justification of Scaling Fermat Distances with Dimension
We propose scale the metric using temperature equal to the dimension of the data. In order justify this, consider the
behavior of a generalized Gaussian distribution of the form p(z) xx exp(-I|><]* /a). It can be shown that B[[[=][]] D and thus the typical distance from the origin of points sampled from p scales as D1/a, The score in this case is s(x) -||0-||-22.
Due to the rotational symmetry of the distribution, we can suppose without loss generality that a sample x lies in the plane of the first two coordinates and hence the problem reduces to the 2-dimensional one. However, since the typical distance of sampled points scales with the dimension, in order to avoid very exaggerated geodesics in high dimensions which deviate sharply towards the origin, we need to scale the density appropriately. By rescaling in Equation (3) by D1/a, the score scales like 1)//a and so the terms in the equation scale as:
h - -3(s 7) + Bs|1111 + DD/aa- 11/0 (.-5) + 3D(a+1)/ 3D(at1)/asl11// s|1>11
(32)
D2 23D(s 4)9 + 3D=
(33)
Therefore the geodesic equation can scale consistently if, and only if, 3 = 1/D. The same reasoning applies to the reparameterized Equation (4).
We can apply a similar argument to the Student-t distribution of the form p(x) x (1 + (1///))-00+ with v > 2. In the limit // + oo. the distribution becomes normal, and otherwise characterized by heavier tails. The expected squared norm ||||=|||| = DD/(v 2). so the typical distance from the origin scales like vD, in a Gaussian distribution. The score function has the form > |=||="; and with x vD, we have x vD, again the same behavior as the Gaussian.
=+D
Hence the scaling behavior is approximately the same as the Gaussian distribution, so the above argument for 3 = 1/D applics as well. Given that the same scaling behavior holds as in the generalized Gaussian, it suggests that 3 = 1/D is natural default choice across a wide range of distributions.
12

---

Learning Distances from Data with Normalizing Flows and Score Matching
Normal
GMM
Circle
Spiral
Spinass
Figure 8. Samples from the two-dimensional datasets.
B. Experimental Details
B.1. Software Libraries
We used the numpy (Harris et al., 2020) and scipy (Virtanen et al., 2020) Python libraries for basic numerical and graph operations and matplotlib (Hunter, 2007) for plotting. We used s scikit-learn (Pedregosa et al., 2011) to compute nearest neighbor graphs and pytorch (Paszke al., 2019) for training neural networks. We used FrEIA (Ardizzone et al., 2018-2022) for constructing normalizing flows.
B.2, Two- Dineensinninnl Datasets
We use the following two-dimensional datasets. All are implemented as Gaussian mixture models. The first two are explicitly Gaussian, the final three are GMMs with 50 components fitted using scikit mesiitllears to data generated by adding noise certain geometric structures.
1. Standard normal distribution
2. Gaussian mixture with 3 components:
Means: #4 = (2, 1.4), #12 = (6.5, 6.3), 43 (8,1)
2.5
33
0
2
-0.8)
Covariances: Σj = 2.5
Σn
Eg
3
3
-0.8
2
Weights: (0.25, 0.5, 0.25)
3.
Circle with radius and Gaussian noise (a = 0.08)
4. Single spiral with 1.75 rotations (3.57 radians) and Gaussian noise (a = 0.05)
5. Two interleaved spirals with rotation each and Gaussian noise (a 0.045)
The spirals reach radius of after single rotation.
13

---

Learning Distances from Data with Normalizing Flows and Score Matching
B.3. Normalizing Flow Training
We used the FrEIA library (Ardizzone al., 2018-2022) design normalizing flows using spline coupling blocks (Durkan al., 2019). We used the following hyperparameters:
Blocks: for two-dimensional data, for standard normal D dimensions
Spline bins: 10
Domain clamping in the splines:
ActNorm between blocks
Noise 10-3 added to the training data
LR x 10-1
Weight decay = 10-6
Fully-connected coupling block subnets with hidden layers of 64 dimensions each, with ReLU activation and BatchNorm
Training iterations (as function of training set size n): 2500 22//1000
Batch size 256 if n 400, else |n/3]
Adam optimizer
B.4. Score Model Training
We use fully connected networks to predict the score, with 6 hidden layers of varying width depending on the input dimension. The architecture consists of:
Input layer - 6 hidden layers - Output layer (dimension matching input)
Hidden layer widths:
3 < D < 128
- 6 < < 170
9 < D < 25: 200
LR = 10-3
Weight decay= = 10-6
Softplus activation function
Training iterations (as function of training set size n): 2500 x 2n//1000
Batch size = 256 if n > 400, else [n/3]
Adam optimizer
In all cases, the chosen model checkpoint was that which achieved the lowest training loss.
B.5. Compute Resources
For the two-dimensional datasets, we performed all experiments on machine with NVIDIA GeForce RTX 2070 GPU and 32 GB of RAM. It took approximately 24 hours to run all experiments.
For the higher dimensional standard normal distributions, we performed all experiments machine with two GPUs: a NVIDIA GeForce RTX 2070 and a NVIDIA GeForce RTX 2080 Ti Rev. A. The machine has 32 GB of RAM. It took approximately hours to run all experiments.
In addition we did some preliminary experiments on the first machine, totalling approximately 48 hours of compute time.
14

---

Learning Distances from Data with Normalizing Flows and Score Matching
C. Additional Experimental Results
C.I. MNIST experiments
Table 1. For each digit, we compute the mean and standard deviation of the log distances between the cluster mean and 20 randomly sampled points that class (intra-class log distance), as well between the cluster mean 20 randomly sampled points other classes (inter-class log distance). We also compute the difference, which gives an indication of how well-separated the digits are in terms log distance.
Digit
Intra-Class Log Distance
Inter-Class Log Distance
Difference (Inter-Intra)
2.41 = 0.40
3.73 # 0..7
1.32 # 0.49
1.69 = 0.49
3.73 = 0.31
2.04 + 0.58
3.56 0.29
4.05 + 0.19
0.49 # 0.35
3.03 0.33
3.78 0.24
0.75 # 0.41
4
3.04 0.38
3.77 0.36
0.73 + 0.53
5
3.43 0.35
4.02 0.18
0.59 = 0.39
6
2.89 # 0.41
3.79 # 0.31
0.91 # 0.51
2.86 # 0.47
3,81 + 0.30
0.95 # 0.56
3.29 = 0.33
3.74 + 0.21
0.46 + 0.39
2.76 + 0.45
3.56 0.38
0.80 + 0.59
3
3
3
3
2
7
66
6
6
L
L
L
2
2
2
5
5
5
5
9
9
9
9
9
9
4
/|
I
1
1
7
7
7
7
77
L
L
2
2
2
2
2
2
2
2
2
Figure 9. Each row visualizes an MNIST geodesic, obtained by relaxation (Algorithm 1) using a learned score model within 10-
dimensional autoencoder latent space.
We train autoencoder to compress MNIST (LeCun et al., 2002) digits to 10-dimensional space and subsequently train sliced score matching model (Song et al., 2020) and normalizing flow (Dinh et al., 2017) on the latent distribution. Both encoder and decoder use fully connected layers with ReLU activations, with a single hidden layer of width 400. The score model has input and output dimension 10 and has hidden layers of width 128, 256, and 128, with SiLU activations. The normalizing flow is built with the FrEIA framework (Ardizzone et al., 2018-2022), made up affine coupling blocks, with permutation between each block. Each block uses coupling network with hidden layers of width 128 and SiLU activations. All models trained for 400 epochs with the Adam optimizer, the autoencoder with learning rate 10-3 and the score model and normalizing flow with learning rate 14
To compute distances between two points sampled from the latent distribution, we the same procedure as in the main text: i) construct a censity-weighter graph (Algorithm 2) from the normalizing flow model ii) use Dijkstra's algorithm to find the shortest path between the points according to the graph iii) solve the relaxation problem (Algorithm 1) with the score model, using the graph-based path as initialization iv) evaluate the length of the relaxed path using the normalizing flow model,
15

---

Learning Distances from Data with Normalizing Flows and Score Matching
according to Equation (2).
In Table I we report mean and standard deviation of log distances between each digit cluster mean and 20 random samples of the same or other digits. We see that same-class distances are always lower than other-class distances, up to a safe margin of error. We also make the following observations: Is cluster tightly, likely due to having fewer active pixels and therefore having higher likelihoods; 8s are more ambiguous, likely due to easy deformations into other digits (it is quite easy to remove some pixels to turn an 8 into a 3, 5, 6 or 9). We leave deeper analysis of Fermat distances on MNIST to future work.
In Figure 9 we visualize geodesics between 5 randomly sampled pairs from the latent distribution. The geodesics are solved in the latent space and the resulting latent points are decoded to image space for visualization. We see that the geodesics follow high-density regions of the distribution, leading to interpretable and reasonable transitions between the digits.
C.2. Convergence of Alternative Graph Edge Weightings
0.
Trae Deusity
Tue Dealty
Men Tog Pail Balio
655
Emvenner of Moan (Tne))
Inverve Mau (Eeiiaated)
Vonu Log Phill Ballo
0.5
Mau Lovee (Beinatol))
Meen hoons (Trns)
64
6.5
0.3-
62-
002
6.1-
0.1
00
00
30P
165
10'
102
1PP
102
Sumplex
Surples
True Demety Mox (Dittmated)
Time Desetr
66-
Mnn (Corinated)
Moan Pail haae
Max (Trac)
Mon twit Rallo
Min (Thue)
6.5
1244
64
63-
0.2
6.-
60
10
16
20'
10
II
10
107
Samples
Sumples
Figure 10. Alteraaive edge weightings based on the nearest neighbor density estimator converge poorly. The same weightings based on the ground truth density converge well. Experiments performed on the GMM two-dimensional dataset.
We try out other edge weights, as described in the main text. Given the nearest-neighhor density estimates p(Xi) o (minm [|X,- - II= for data points Xt, we approximate the graph edges by Gtm = ||X, - xn ( + x2m) where pis density estimate for the midpoint. We use 4 variants:
1. Inverse of Mean: p((X) + Xm)) XXm))) = (f(x)) + ((X
2. Mean of Inverse: 1//((XI + xm.)/2 = (1/((X)) + 1//(2,)))
3. Max: p((X,+ + km)/2)= = max(()()),,
4. Min: P((X, + X2) m(i()()))..(Xmm)))
We find that all have very similar performance to the power-weighted graph, whereas the equivalent estimators using the ground truth density quickly converge. See Figure 10.
C.3. Trfde-off Between Log Density and Score Estimation
We speculate that there is a fundamental trade-off between accurately estimating log densities and their gradients (score functions) that applies to all density estimation methods, including the normalizing flows used in our main experiments. To
16

---

Learning Distances from Data with Normalizing Flows and Score Matching
Density
Error
123
Imgraae Squard
Figure 11. Mean integrated squared error (MISE) for log density and score estimation using KDE with varying bandwidth. The optimal bandwidth for log density estimation is smaller than the optimal bandwidth for estimation.
illustrate this principle in a simple setting where we can compute exact errors, we analyze kernel density estimation (KDE) of the standard normal distribution in one dimension using Gaussian kernels, with 1000 samples drawn from the distribution. While this example is much simpler than the high-dimensional problems and neural network models considered in the main text, it clearly demonstrates how optimizing for density estimation can lead to poor gradient estimates and vice versa.
For set of samples {z,) the KDE estimate with bandwidth is:
exp (- (x- - 22,)
p(x) no 2
(34)
2n²
The log density and score estimates are then:
(x - x;
log p(x) (og(no//2)) + log Σ
exp
(35)
202
i=1
En (1r x1:) exp
2
8(x)
(36)
a² E exp
(x-2,)² 20
As shown in Figure 11, the optimal bandwidth for minimizing the mean integrated squared error (MISE) of the log density estimate is smaller than the optimal bandwidth for minimizing the MISE of the score estimate. This reflects a fundamental trade-off: accurate density estimation requires capturing fine details of the distribution, while accurate score estimation requires smoothing out noise to obtain reliable gradients.
Figure 12 visualizes this trade-off. When using the bandwidth optimal for log density estimation, the score estimates are quite noisy. Conversely, when using the bandwidth optimal for score estimation, the log density estimates are oversmoothed. This observation helps explain why we obtain better results using dedicated score models rather than differentiating density models.
17

---

Learning Distances from Data with Normalizing Flows and Score Matching
Log Density Comparison (0=D.278) Optimal Log Density
Log Density Comparisan (o-0.359) Optimal
Log Donsity
Log Dersity
-8
KDE log censtyy Tue log densty
KDE log density
Tue to denstty
3
2
Score Comparison (a=0.278) Optimal for Log Density
Score Comparison (p=0 0.359)
Optimal for Score
KDE
Score
Score
-10
-4
-2
Figure 12. Comparison KDE estimates using different bandwidths. Top row: Log density estimates. Bottom row: Score estimates. Left column: Using the bandwidth that minimizes log density MISE. Right column: Using the bandwidth that minimizes score MISE. Ground truth shown in dotted orange, KDE estimate solid blue. The optimal bandwidth for log density leads to noisy score estimates, while the optimal bandwidth for estimation leads oversmoothed log densities.
18


Reviews:
[{"review_id": "VHfhJcu84n", "reviewer": "Reviewer_AqTp", "summary": "In this paper, the authors propose to learn a Riemannian metric from data using a class of Fermat metrics, which are metrics that are equal to the Euclidean metric rescaled at each point by (a power of) the reciprocal of the probability density of the data. This way, geodesics tend to follow high density regions, which is desirable e.g. to interpolate between data points. The authors point out limitations of current estimators of this type of distances, mostly because of poor density estimation and the curse of dimensionality. They introduce a relaxation scheme to numerically solve for geodesic segments linking two fixed endpoints when the ground truth density is known (this allows computing geodesics for known densities). They also modify existing graph based methods by estimating the density with normalizing flows. They finally introduce another method which combines the relaxation scheme with score estimation. The methods are tested on 2D examples whose ground truth densities are known, and the score based method is tested on higher dimensional standard Gaussians. The competing methods are the graph based methods that don’t use deep learning techniques for density estimation.\n\n## Update after rebuttal\n\nThanks to the authors for the responses and the preliminary experiments on MNIST. I would indeed be nice to add an appendix to discuss these, though more extensive experiments and theoretical results would have strengthened the paper even more. I still have an overall positive opinion on this paper, so I will maintain my score.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Rebuttal to Reviewer AqTp\n\nWe thank the reviewer for their thoughtful and constructive feedback. Below, we address each of the main points and questions raised.\n\n---\n\n### How large can the dimension become for the relaxation method to stay efficient and tractable? The finite difference scheme must require finer discretization in higher dimensions, I suppose?\n\nThe tractability of the relaxation method primarily depends on the complexity of the data, rather than the ambient dimension. For example:\n\n- Uniform distributions: Our method correctly converges to straight lines even with very coarse discretization.\n- Standard normal: Reduces to the 2D case (due to rotational symmetry when using β = 1/D), and we find that 20 segments are sufficient for accurate geodesics, even in high dimensions.\n\nMore complex distributions may require finer discretization for accuracy, but there is no intrinsic limitation due to dimension itself. In contrast, the NF graph method becomes increasingly ineffective in higher dimensions, as filling the space requires exponentially more samples due to the curse of dimensionality.\n\n---\n\n### How would the method perform in small and simple image datasets such as MNIST or CIFAR?\n\nWe agree that applying the method to image data is an exciting next step. While ground truth densities are unavailable, qualitative insights (e.g., visualizing shortest paths) can still be informative.\n\nAs a first step, we train an autoencoder to compress MNIST to a 10-dimensional latent space, then train a normalizing flow and score model on the latent distribution. Solving for geodesics yields plausible digit interpolations.\n\nAs an initial quantitative measure, we report mean and standard deviation of log distances between each digit cluster mean and 200 random samples of the same or other digits:\n\n| Digit |     Same Class     |    Other Class     | Diff (Other-Same)  |\n|-------|--------------------|--------------------|--------------------|\n|   0   |    2.41 ±  0.40    |    3.73 ±  0.27    |    1.32 ±  0.49    |\n|   1   |    1.69 ±  0.49    |    3.73 ±  0.31    |    2.04 ±  0.58    |\n|   2   |    3.56 ±  0.29    |    4.05 ±  0.19    |    0.49 ±  0.35    |\n|   3   |    3.03 ±  0.33    |    3.78 ±  0.24    |    0.75 ±  0.41    |\n|   4   |    3.04 ±  0.38    |    3.77 ±  0.36    |    0.73 ±  0.53    |\n|   5   |    3.43 ±  0.35    |    4.02 ±  0.18    |    0.59 ±  0.39    |\n|   6   |    2.89 ±  0.41    |    3.79 ±  0.31    |    0.91 ±  0.51    |\n|   7   |    2.86 ±  0.47    |    3.81 ±  0.30    |    0.95 ±  0.56    |\n|   8   |    3.29 ±  0.33    |    3.74 ±  0.21    |    0.46 ±  0.39    |\n|   9   |    2.76 ±  0.45    |    3.56 ±  0.38    |    0.80 ±  0.59    |\n\nWe see that same-class distances are always lower than other-class distances, up to a safe margin of error. We also note: 1s cluster tightly, likely due to fewer active pixels and therefore higher likelihoods; 8s are more ambiguous, likely due to easy deformations into other digits (it is quite easy to remove some pixels to turn an 8 into a 3, 5, 6 or 9). We will include these results and visualizations in the paper, but leave deeper analysis of Fermat distances on MNIST to future work.\n\n---\n\n### I assume the geodesic equation for Fermat metrics is a known result (reference?), and that only the reparameterization is new.\n\nYes, the form of the geodesic equation for conformal metrics follows from standard Riemannian geometry (e.g., see Appendix G of *Spacetime and Geometry*, Carroll, 2004). The reparameterization to constant Euclidean speed is indeed new and key to our numerically stable relaxation scheme.\n\n---\n\n### I suggest adding a plot showing an example of geodesics recovered on a dataset by the different methods.\n\nAgreed—we will add such a figure to the main text to complement Figure 4, comparing several methods on the same dataset.\n\n---\n\n### A plot with approximated geodesics with various discretization levels would be interesting.\n\nWe will add such a figure to the appendix to illustrate convergence behavior and practical trade-offs.\n\n---\n\n### Minor Comments\n\n- Figures 1 and 4: Will ensure they are clearly referenced.\n- Repeated sentence on p3: Will be removed.\n- Fig. 10: Will standardize y-axis ranges for better visual comparison.\n\n---\n\nWe thank the reviewer again for their detailed review and insightful suggestions. We believe these improvements will significantly strengthen the clarity and impact of the paper.", "claims_and_evidence": "The paper makes a certain number of claims compared to the existing literature. First the capacity to use the relaxation algorithm for datasets with known densities enables quantitative comparisons, which were not extensively performed previously. This in particular can help show shortcomings of classical graph based methods in the experiments. The density estimation for graph based methods seems to be the issue, and indeed using NF for this purpose improves performance. The combination of the score estimation and the relaxation method seems to scale better with dimensionality(up to d = 30), albeit only a single experiment using standard normals is performed. \nOf course, the claim of the score based method to scale well with dimensionality calls for experiments on e.g. image datasets. Score based methods are known to scale very well with dimensionality, so this raises two questions : \n1) How large can the dimension become for the relaxation method to stay efficient and tractable? The finite difference scheme must require finer discretization in higher dimensions, I suppose? Maybe in that case the NF graph based method would be a better choice if there are sufficiently many samples?\n2) How would the method perform in small and simple image datasets such as MNIST or CIFAR (possibly working in a latent space if the input dimension of the images turn out to be too large) ? I understand that in this case no GT density is available but qualitative behavior can be analyzed and it would be interesting to see shortest paths in image space.", "methods_and_evaluation_criteria": "Though I am not very familiar with the literature on learning Fermat metrics, the introduction of a method that allows to compute ‘GT’ geodesics for known densities is definitely a good thing for the domain if it was hard to compare to true geodesics before. It apparently helped identifying shortcomings of previous classical methods, in spite of apparently reassuring theoretical guarantees.", "theoretical_claims": "I checked the proofs and derivations. I assume the geodesic equation for Fermat metrics is a known result (reference?) and that only the reparameterization is new. This should be precised in the text. \nIs there an existing result on the convergence of the relaxation scheme (when successive values become closer and closer) ? This seems like an important result to have to further strengthen the method. Similarly, investigating further the reason why convergence of the classical graph based method is so low in practice would be interesting.", "experimental_designs_or_analyses": "The experiments are sound and show the improvements brought by the two methods that introduced. The remarks on the difficulty of estimating the density and its score simultaneously are interesting and are somewhat reminiscent of Heisenberg’s uncertainty principle.", "relation_to_broader_scientific_literature": "I am not an expert on Fermat metric learning, so I cannot comment on the exhaustivity of the related work on methods for this. However, the paper is clearly positioned among the papers that are cited.", "essential_references_not_discussed": "No missing references to my knowledge.", "other_strengths_and_weaknesses": "- The lack of experiments in higher dimensions is a bit of a shame, as it would really tell about the potential of the method for generic ML tasks and datasets.\n- I suggest adding a plot showing an example of geodesics recovered on a dataset by the different methods. Something similar to Figure 4, but not restricted to graph based methods.\n- In line with my comment on the convergence of the relaxation scheme, a plot with approximated geodesics with various discretization levels would be interesting (how is this chosen in practice?).\n\nOverall the paper is interesting with a number of compelling elements, but it feels slightly limited in terms of theoretical guarantees and experiments on higher dimensional datasets would be most welcome.", "other_comments_or_suggestions": "- Figs 1 and 4 are not referenced in the text.\n- There is a repeated sentence in p3, right column, below Fig. 2.\n- I suggest setting the same extent on the y axes of both plots for each columns of Fig. 10 for better comparison between both cases.", "questions_for_authors": "See my questions the boxes above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "3Oe4dBqxzk", "reviewer": "Reviewer_EAyN", "summary": "This paper addresses the problem of learning distance metrics from data, specifically focusing on density-based distances (DBDs). The authors highlight that existing methods for estimating Fermat distances suffer from poor convergence and scaling issues in high dimensions due to inaccurate density estimates and insufficiently smooth geodesics. To tackle these challenges, the paper introduces two main improvements: using normalizing flows to learn more accurate densities and refining geodesics with a learned score model. Additionally, the authors propose a dimension-adapted Fermat distance to improve scaling and numerical stability in higher dimensions. The core idea is to leverage density/score estimation using modern deep learning to improve the practical applicability of density-based distances.\n\n## update after rebuttal\nMy questions and concerns have been addressed and I would be very happy to see the paper accepted.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Rebuttal to Reviewer EAyN\n\nWe thank the reviewer for their detailed and constructive comments.\n\n---\n\n### Simplicity of Experiments\n\nPlease see our response to reviewer AqTp for results on an experiment on MNIST.\n\n---\n\n### Missing Recent Work\n\nThank you for pointing out our oversight. We will add both references to the discussion of related work.\n\n---\n\n### Time Complexity and Scaling\n\nScaling the dimension does not necessarily increase the number of segments needed in the relaxation algorithm — this depends on the complexity of the dataset (see also our response to Reviewer AqTp). The complexity per iteration of the loop scales linearly with dimension, due to norm and inner product computations. So overall, the complexity of relaxation is roughly linear in dimension.\n\nFor the simple case of a uniform distribution (where the score is zero and the solution is a straight line), we can show that convergence requires $O(n^2)$ iterations, where $n$ is the number of segments. While a full analysis in general cases is more challenging, we can provide this estimate as a starting point and are happy to add more details in the appendix if helpful.\n\nThe graph-based algorithms also scale linearly in dimension, but the number of nearest neighbors may need to increase with dimension to preserve local structure. However, as shown in Figure 6, graph-based methods degrade in high dimensions due to sparsity and noise, often before computational scaling becomes the limiting factor.\n\n---\n\n### Fermat Distance Scaling in Non-Gaussian Distributions\n\nYou are correct — our original discussion was overly focused on the Gaussian case. For a generalized Gaussian distribution of the form $p(x) \\propto \\exp(-\\|x\\|^\\alpha / \\alpha)$, it can be shown that $\\mathbb{E}[\\|x\\|^\\alpha] = D$, and thus the typical distance from the origin scales as $D^{1/\\alpha}$. The score in this case is $s(x) = -\\|x\\|^{\\alpha - 2} x$.\n\nBy rescaling $\\gamma$ in Eq. (3) by $D^{1/\\alpha}$, all terms in the equation scale consistently if and only if $\\beta = 1/D$. The same reasoning applies to the reparameterized equation (Eq. (4)).\n\nA similar argument can be made for the Student-t distribution, where $p(x) \\propto \\left(1 + \\|x\\|^2 / \\nu\\right)^{-(\\nu + D)/2}$. The scaling is approximate in this case but still suggests that $\\beta = 1/D$ is a natural default across a wide range of distributions.\n\nThank you for encouraging us to clarify and strengthen this argument — we will update the text and appendix accordingly.\n\n---\n\n### Speed and Initialization of Algorithm 1\n\nPlease see our earlier comments on time complexity. Regarding initialization: it is indeed possible to initialize the relaxation algorithm using simpler paths, such as power-weighted shortest paths based only on Euclidean distances. This often works well in practice and is especially useful when only the score is available (without a density model).\n\n---\n\n### Computation of Ground Truth Geodesic Distance (Eq. 13)\n\nWe use Algorithm 1 with the ground truth score function to compute true geodesics. Initial trajectories are obtained using a graph-based approximation. This is discussed in Section 3.1.\n\n---\n\n### Absence of Graphs for NF-Based Relaxation\n\nYou are right — we do not include plots of NF-based relaxation because the scores derived from differentiating normalizing flow models are quite noisy. This frequently causes the relaxation algorithm to diverge, making the results unreliable for evaluation. This is precisely why we trained separate score models.\n\n---\n\n### Hybrid Training of NFs with Score Matching\n\nWe experimented with combining maximum likelihood training and sliced score matching. Unfortunately, this either did not improve the quality of the scores or led to training instability. We are not certain of the root cause, but we suspect it reflects the trade-off between optimizing for density accuracy and gradient (score) accuracy — see the discussion in Appendix C.2.\n\n---\n\n### Line 434: Unifying Normalizing Flows and Score Models\n\nBy this, we mean developing a single model that yields both accurate densities and scores. Our attempts to do so with normalizing flows were not successful (see comment above), so we leave this as a challenge for future work.\n\n---\n\n### Line 102: Numerically Stable Relaxation Method\n\nWe refer here to our relaxation method based on a reparameterized geodesic equation, which enforces constant Euclidean speed. This improves numerical stability and simplifies discretization. To our knowledge, this reparameterization has not appeared in prior work and is a novel contribution of this paper.\n\n---\n\nWe thank the reviewer again for their thoughtful feedback and for helping us sharpen the clarity and scope of the work.", "claims_and_evidence": "The claims made in the submission are generally well-supported by clear and convincing evidence. The authors provide both theoretical grounding and empirical results to back up their contributions.  \n* One key claim is that previous graph-based methods exhibit poor convergence.  This is supported by the experimental results in Section 4.1.1 and Figure 5, which show that the power-weighted graph method converges slowly.   \n* The claim that normalizing flows improve density estimation is supported by the improved convergence rates observed when using normalizing flows for edge weights (Section 4.1.3 and Figure 5).  \n* The effectiveness of score-based relaxation in higher dimensions is demonstrated in Section 4.2 and Figure 6, where the relaxation method maintains performance while graph-based methods degrade.   \n\nHowever, there still remains a large gap as most of the examples considered in this paper are of extremely simplistic distributions and lower dimensions (compared to realistic datasets used in most fields).", "methods_and_evaluation_criteria": "The datasets used appear rather simplistic and hand-crafted. For example, none of the datasets appear to have disconnected support.\nLPR appears to be a good performance metric, but following other similar works, I would suggest also considering geodesics standard geodesic error and geodesic variation errors as in https://arxiv.org/pdf/2403.06612 (eqns 115, 116).", "theoretical_claims": "I have checked the correctness of the proofs for the theoretical claims in appendix A1. The derivation of the geodesic equation in Theorem A.2 and the reparameterization in Theorem A.3 appear to be correct, with steps clearly laid out.", "experimental_designs_or_analyses": "They appear sound.", "relation_to_broader_scientific_literature": "The paper appears to miss the recent proposals for learning geometry from data using score models (https://arxiv.org/abs/2405.14780) and normalizing flows (https://arxiv.org/abs/2410.01950).", "essential_references_not_discussed": "see above.", "other_strengths_and_weaknesses": "* The paper in its current form does not seem to address the question of time complexity required - especially how it relates to the dimension scaling. \n*", "other_comments_or_suggestions": "* I would strongly suggest clearly separating which methods are graph based and which are not, as currently it is not clear at times.\n\nand see below.", "questions_for_authors": "* Fermat Distance dimension scaling - the example provided makes a lot of sense when considering gaussian distributions. But it is unclear to me whether the same scaling should be used when considering a different distribution - could you comment/elaborate on validity in other cases? \n* Line 133 - you mention that beta = 1 is common in previous works, which would be useful to have a couple of references for. \n* Could you comment on how slow Algorithm 1 is? Also, if the initialization is not done using normalizing flows, is the geodesic still computed accurately?\n* How is the true geodesics distance in (13) computed?\n* Line 358 - when referring to supplementary materials, please provide where in supplemetaries this can be found. \n* There appear to be no graphs for using normalizing flows directly for relaxation. Could you either point me to where they can be found or include them? You appear to claim that they are not very good and it would be good to understand what exactly you mean by that. \n* Building on the previous question - NF trianed with maximizing likelihood would make sense to not be very good at approximating the score. Have you tried training an NF with a mix of likelihood and sliced score matching? It is not clear to me why both can be optimized for.\n* Line 434 - what do you mean \"Unifying nf and score models into a single framework?\"\n* Line 102 - you claim to \"Introduce a numerically stable relaxation method\" - what exactly do you mean by this? Is this your novel contribution?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "m6861FBqBX", "reviewer": "Reviewer_XLew", "summary": "The paper presents a method to learn distances from data by integrating normalizing flows and score matching into the computation of density-based distances (DBDs), specifically Fermat distances. It addresses the shortcomings of existing methods by introducing a stable numerical approach to compute true geodesics through normalizing flows and enhancing the smoothness of trajectories in high dimensions using score matching. The authors validate their approach by demonstrating faster convergence and improved accuracy over traditional graph-based methods, particularly in complex Gaussian mixtures and higher-dimensional spaces.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Rebuttal to Reviewer XLew\n\nThank you for your positive assessment of our work and for your thoughtful suggestions.\n\n---\n\n### Preliminary Results on Real-World Dataset\n\nPlease see our response to reviewer AqTp for a preliminary experiment on MNIST. We find that the distances obtained with our method behave as expected, and yield interesting insights into the relationships between the different classes.\n\n---\n\n### Computational Efficiency Compared to Graph-Based Methods\n\nPlease see the detailed discussion in our response to Reviewer EAyN. In short, the relaxation method (Algorithm 1) has per-iteration complexity that scales linearly with dimension due to the need to compute norms and inner products. The number of segments needed depends more on the complexity of the data than on the ambient dimension.\n\nGraph-based methods similarly scale linearly with dimension in terms of distance computations, but their performance degrades significantly in higher dimensions, as shown in Figure 6, due to the sparsity of samples and the limitations of nearest-neighbor approximations.\n\n---\n\n### Assumptions Underlying Theoretical Results\n\nThe theorems provided in the paper, particularly the derivation of the geodesic equations and their reparameterization, rely on standard smoothness assumptions:\n\n- The density function $p(x)$ is assumed to be differentiable and non-zero in the region of interest.\n- The score function $s(x) = \\nabla \\log p(x)$ is assumed to be Lipschitz continuous.\n- The conformal metric defined by $p(x)$ is smooth.\n\nThese conditions ensure the existence of solutions to the geodesic equations and support the stability of the relaxation method. They are commonly satisfied in practice, particularly when using neural networks for score estimation.\n\n---\n\nThank you again for your encouraging review and support.", "claims_and_evidence": "The claims in the paper are supported by a mix of theoretical development and experimental results. The authors present a clear improvement in the accuracy and computational feasibility of estimating Fermat distances using their method. However, the claim regarding the general applicability of their approach in \"real-world, high-dimensional applications\" could be seen as overly optimistic, given the experiments are primarily on synthetic or simplified datasets. More evidence from diverse and complex real-world datasets would strengthen this claim.", "methods_and_evaluation_criteria": "The methods proposed, including the use of normalizing flows for accurate density estimation and score matching for refining geodesics, are well-suited to the stated problem of improving the computation of DBDs in high-dimensional spaces. The evaluation criteria, particularly the use of Log Path Ratio (LPR) for comparing different methods, is relevant and provides a clear metric for assessing improvements over existing approaches.", "theoretical_claims": "The paper outlines several theorems related to the geodesic equations and their solutions using the proposed methods. The derivation of these theorems, such as the constant Euclidean speed parameterization and its impact on the stability of numerical methods, appears logically sound. No specific issues were identified in the proofs provided, but a more detailed external validation of these theoretical claims would ensure their robustness.", "experimental_designs_or_analyses": "The experimental design seems valid for demonstrating the advantages of the proposed methods over traditional approaches. The use of synthetic datasets, while controlled, is appropriate for illustrating the performance improvements in known environments. The methodological setup, including comparisons to ground truth geodesics and other baseline methods, allows for a clear demonstration of the proposed method's superiority in terms of convergence rates and accuracy. However, expanding the experiments to include a broader range of real-world datasets would help validate the practical applicability of the methods outside controlled experimental conditions.", "relation_to_broader_scientific_literature": "The key contributions of the paper, specifically the integration of normalizing flows and score matching into the computation of density-based distances, draw significantly from the established fields of geometric deep learning and statistical machine learning. Prior research in normalizing flows has typically focused on improving the accuracy and efficiency of probability density function estimation in complex data distributions (e.g., Rezende and Mohamed, 2015). The application of these flows to compute Fermat distances introduces a novel intersection of geometric learning with probabilistic modeling, expanding upon works like those by Arjovsky et al. (2017) in Wasserstein GANs that explore distance metrics in latent spaces. The use of score matching, introduced by Hyvärinen (2005), for smoothing trajectory calculations in high-dimensional spaces further builds on the idea of refining probabilistic estimations without explicit density estimations.", "essential_references_not_discussed": "Beneficial not essential \n\nWhile the paper adequately cites foundational works on normalizing flows and score matching, it may lack references to some pertinent studies that bridge these concepts more directly with geometric applications in machine learning. For instance, research by Bronstein et al. (2017) on geometric deep learning frameworks could provide additional context for the application of geometric principles in learning tasks. Additionally, recent advancements in computational geometry for machine learning, such as those presented in the works by Memoli (2011) on Gromov-Wasserstein distances, could enhance the theoretical underpinnings and practical applications of the methods discussed in this work.", "other_strengths_and_weaknesses": "Strengths:\nOriginality: The paper's approach to integrating normalizing flows with score matching to compute Fermat distances is highly original. This creative combination of techniques from different areas of machine learning could set a new precedent in the field.\nSignificance: The potential impact of this method in improving the computational feasibility and accuracy of distance calculations in high-dimensional spaces is significant, especially for applications in complex data analysis and geometric learning.\nClarity: The paper is well-written, with clear explanations of the methods and their theoretical foundations, making complex concepts accessible to readers.\n\nWeaknesses:\nGeneralizability: The paper primarily demonstrates results on synthetic datasets. The generalizability of the approach to real-world datasets and its performance in truly unstructured environments remain to be validated.\nComplexity Discussion: The computational complexity and resource demands of the proposed methods are not thoroughly discussed, which could be crucial for practical applications needing scalability considerations.", "other_comments_or_suggestions": "I really liked Figure 1; it offers immediate clarity and a high level overview of your work.\n\nWhile reviewing the paper, I noticed a few minor typographical errors that could be corrected to enhance the overall readability and professionalism of the manuscript:\n\nOn page 3, second paragraph, \"theorm\" should be corrected to \"theorem.\"\nOn page 5, in Figure 2's caption, \"illustraiting\" should be changed to \"illustrating.\"\n\nIt would also be beneficial to include a supplementary section with additional details on the parameter settings for the normalizing flows and score matching algorithms used in the experiments. This addition would aid in replicating the results and understanding the sensitivity of the proposed method to different configurations.", "questions_for_authors": "1. Given that the paper primarily focuses on synthetic datasets, can you provide insights or preliminary results on how the proposed method performs on real-world datasets, particularly those with noise and irregular distributions?\n\n2. Could you elaborate on the computational efficiency of your method, especially in comparison to traditional graph-based methods? What are the traditional graph-based methods that you may choose as a baseline? Specifically, what are the computational costs in terms of time and resources when applied to larger datasets?\n\n3. Are there any specific assumptions or conditions under which the proposed theorems hold? How do these assumptions affect the generalizability of your findings?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ntSEPOveIc", "reviewer": "Reviewer_nn54", "summary": "The article proposes and compares several methods for estimating distances derived from Riemannian metrics that reflect the data distribution. In particular, the chosen metric should \"compress\" distances in regions of high mass concentration and \"stretch\" distances where the mass is lower. To achieve this, the authors consider conformal metrics (i.e., proportional to the identity) that are inversely proportional to the density of the underlying distribution of the observations. At a point $x$, the density-based metric is given by $\\frac{I_d}{p(x)^{\\frac{2}{d}}}$ where $p$ is the density of the observations and $d$ is the dimension of the space. Given a dataset where the density is unknown, the challenge is to estimate density-based distances and geodesics.\n\nThe authors propose two methods for estimating a density-based geodesic between any two points:\n\n1) The first method relies on numerically solving (through relaxation) the differential equation whose solutions are geodesics. The authors provide an expression for this equation in the case of the density-based metric, which explicitly involves the derivative of the density. There are two approaches to estimate this quantity:\n      - Computing the derivative of the density numerically, where the density itself is estimated from the data using a normalizing flow.\n      - Directly estimating the density derivative using score matching—a method ultimately chosen because it is more efficient.\n\n2) The second method approximates a density-based geodesic by finding the shortest path in a graph where edge weights are proportional to the density-based distance. This distance can be determined in two ways:\n      - Using the method of Bijarl et al. (2012), which estimates the density-based distance as a power of the Euclidean distance. However, the numerical simulations conducted by the authors do not yield good results with this technique.\n      - Estimating the density using a normalizing flow and then plugging it into the density-based metric expression.\n\nFor each of these two methods, the algorithm is detailed in the paper (Although I don't think they are new in the literature). To improve the numerical stability of the relaxation method, geodesics are parametrized to have constant Euclidean speed rather than constant speed relatively to the density-based metric.    \n\n\nThe performance of a method is measured using the log of the ratio between the estimated distance and the true distance (log path ratio). Experiments are conducted on five datasets sampled from probability distributions with known densities (and thus known density-based metrics and distances). Several observations are made :\n- The accuracy of density estimation is crucial for the precision of graph-based method.\n- The relaxation method is much more accurate when the density derivative is estimated directly using score matching.\n- As the dimension increases, the graph-based method suffers from the curse of dimensionality, whereas the relaxation method maintains good performance in high dimensions.\nGiven these observations, the authors recommend using the graph-based method to initialize the relaxation method.\n\n## update after rebuttal\nI thank the authors for their responses and clarification. The modifications announced by the authors seem reasonable and will help clarify the theoretical context of the paper. The topic of the paper, as well as the algorithms proposed by the authors, appear very interesting to me, although the experiments seem insufficient (especially considering that no theoretical guarantees are provided).", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Rebuttal to Reviewer nn54\n\nWe thank the reviewer for their careful and detailed reading of our work, and for the helpful questions and clarifications.\n\n---\n\n### Novelty of the Algorithms\n\n> The algorithms used by the authors are already present in the literature.\n\nWhile the Fermat distance itself is not a novel concept, our algorithms for computing it are, to our knowledge, new. In particular, Algorithm 1 introduces a reparameterization of the geodesic equation that ensures constant Euclidean speed along the curve, making the relaxation tractable and numerically stable. This reparameterization is key to our ability to exactly solve for geodesics in the Fermat metric — a setting where prior work has only used approximations.\n\nAlgorithm 2 also includes novel elements. While graph-based shortest paths are standard, our method departs from previous approaches that rely on Euclidean distances or kernel density estimates to define edge weights. Instead, we use edge weights derived from a learned density model, allowing for more accurate adaptation to complex, high-dimensional data. To our knowledge, this specific use of learned density-driven edge weights in computing Fermat distances has not been explored in the literature.\n\n---\n\n### Metric Tensor Notation\n\n> Does $g$ refer to a metric tensor field ($g: M \\to TM$), which is a continuous function, or does it represent $g = g_p = g(p)$, the inner product defined on $T_p$ for $p \\in M$?\n\nIt is the latter. We will update the notation to make this clear and avoid confusion. Specifically, we will clarify that $g_p$ is a smoothly varying inner product on each tangent space $T_p$.\n\n---\n\n### Scope of the Analysis and Manifold Assumptions\n\n> More generally, I am struggling to understand in which space the analysis in the paper is conducted.\n\nWe apologize for the confusion. You are correct — the analysis is conducted entirely in $\\mathbb{R}^D$ equipped with a conformal Riemannian metric. We do not work on general manifolds. Thank you for pointing out that our use of manifold terminology in Section 2.1 may suggest otherwise. We will revise the exposition to streamline the definitions and make it explicit that our setting is $\\mathbb{R}^D$.\n\n---\n\n### Clarification on “Distance Between Points in a Distribution”\n\n> Does ‘define distance between points in a distribution’ mean ‘define distance between points sampled from a distribution’?\n\nYes, that is correct. Thank you — we will update the phrasing accordingly to avoid ambiguity.\n\n---\n\n### When Are Geodesics Distance-Minimizing?\n\n> Do you know in which Riemannian manifold the geodesics are always distance minimizers?\n\nThis is indeed a challenging problem, and in general there are cases where multiple geodesics exist between two points, some of which are only local minimizers of the path length. For instance, in the case of the circle distribution (see Fig. 7), the shortest geodesic between two nearby points will traverse a small arc, but there is always another valid geodesic that takes the longer path around the circle. If the relaxation algorithm is poorly initialized, it may converge to such a non-minimizing solution. For this reason, we emphasize the importance of using a graph-based shortest path to provide a good initialization. This is described in Section 3.3 (line 289, column 2).\n\n---\n\nWe appreciate the reviewer's in-depth engagement with the technical aspects of the paper and will incorporate the suggested clarifications to improve accessibility and precision.", "claims_and_evidence": "The claims are based on a comparison of methods through experiments. Each performance value is obtained by averaging the log path ratio over 1000 pairs of points selected uniformly at random.\nFor the graph-based method, the performance of the two density estimation methods is compared using estimations on five datasets sampled from five different distributions in $R^2$. For the normal distribution only, a comparison is made between the graph-based method and the relaxation method as a function of the dimension.\nThe authors claim to have highlighted a gap between theory and practice, since their numerical estimates do not achieve the performance guaranteed by the results of Hwang et al. (2016), published in The Annals of Applied Probability. They speculate that the unknown constant in Hwang et al.'s result is responsible for this discrepancy.\nIt would have been interesting to discuss to what extent the assumptions of Hwang et al.'s result are satisfied in the case of the distributions used in the experiments.", "methods_and_evaluation_criteria": "See the paragraph above", "theoretical_claims": "The theoretical claims consist in the development of the geodesic equation for the case of conformal metrics (i.e., proportional to the identity). These equations are provided first with constant speed relative to the density-based metric, and then relative to the Euclidean distance. The proof involves computing the Christoffel symbols and is provided in the appendix.", "experimental_designs_or_analyses": "The experimental methodology is detailed in the appendix. The two deep learning models (normalizing flow and score matching) are trained on a training set (of unknown size), and the retained parameters are those that minimize the loss on the training set.\nThe experimental results are presented as the average of the performance metric over 1000 realizations. No additional indicators are provided (in particular, no empirical standard deviations).", "relation_to_broader_scientific_literature": "The algorithms used by the authors are already present in the literature. The authors suggest using deep learning methods to improve the performance of these algorithms by estimating quantities that are then plugged into the existing methods.\n\nThe approximation of geodesics by the shortest path in a graph, where edge weights depend on the density, was developed by Bijral et al. (2012). Dijkstra’s algorithm (1956) is used to compute the actual shortest path once the edge weights are determined. The authors suggest using a normalizing flow method to compute these weights instead of the approach by Bijral et al. (2012).\n\nThe approximation of geodesics by computing an approximate solution to the geodesic equation is likely already present in the literature, although no citation is provided in the paper. To improve numerical stability in computing an approximate solution to the geodesic equation, the authors seek a solution to the reparametrized equation so that the geodesic has constant Euclidean speed. Additionally, they propose using a score matching method to estimate one of the parameters of the geodesic equation (the gradient of the density).\n\nThe deep learning methods (normalizing flow and score matching) are well-referenced in the article.", "essential_references_not_discussed": "References that allow for the contextualization of the results are cited.", "other_strengths_and_weaknesses": "I am having trouble understanding the given definition of the metric tensor, particularly the sentence: \"The metric tensor g is a smoothly varying positive-definite bilinear form on the tangent space.\" Does $g$ refer to a metric tensor field ($g: \\ M \\rightarrow TM$), which is a continuous function, or does it represent $g = g_p = g(p)$, the inner product defined on $T_p$ for $p \\in M$? In the latter case, it would be more accurate to write $g_p$ and then $∣∣v∣∣_p=g_p(v,v)$ at the beginning of the following paragraph.\n\nMore generally, I am struggling to understand in which space the analysis in the paper is conducted. Section 2.1 suggests that it applies to an arbitrary smooth manifold equipped with a Riemannian metric. If that is the case, what does the probability density $p$  introduced on line 139 correspond to (density with respect to which measure, the Riemannian volume measure?) Similarly, does the change of variable formula (line 255) remain valid for random variables taking values in an arbitrary smooth manifold, where the considered measure is not the Lebesgue measure? Perhaps these objects are commonly used in Riemannian geometry for certain well-characterized spaces, but they should probably be justified more explicitly for readers who are not specialists in the field. \nIn Algorithm 1 (lines 282–284), you sum two elements of the space—what does this mean if the space is, for example, the sphere?\n\nGiven these remarks, perhaps your analysis is only intended to be valid in $R^n$ equipped with a Riemannian metric. If that is the case, I believe the introduction to Riemannian manifolds is unnecessary and could be replaced with the definition of a Riemannian metric on $R^n$ (without having to introduce the notion of tangent spaces).", "other_comments_or_suggestions": "1) L.109 column 2: Does “define distance between points in a distribution” mean “define distance between points sampled from a distribution”?\n2) L.130 column 2: The first paragraph of “Fermat distances which scales with dimension” is written twice.\n3) L.195 column 1: Equation (3) is true for geodesic trajectory so the sentence before is supposed to be “for a geodesic trajectory $\\gamma$ we have”", "questions_for_authors": "Do you know in which Riemannian manifold the geodesics are always distance minimizers? (and so the geodesic equation (Eq. (3), line 197) has as solutions only geodesics that minimize distance). For example, this is not true in the sphere equipped with its usual metric, but perhaps it holds in $R^n$ with a conformal metric (it is true for the identity metric, and maybe conformality preserves this property?).\n\nIf this is not the case, perhaps initializing with the weighted graph method would help the relaxation method converge toward a distance-minimizing geodesic (since the weighted graph method is based on finding the shortest path rather than a on criterion related to the acceleration of the trajectory).", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper deals with the definition and estimation of shortest paths between points in Euclidean space, relative to some unknown distribution. That is, given iid data points in $\mathbb R^d$ with some unknown density $f$, a shortest path between any two points $a,b\in\mathbb R^d$ should avoid regions of low density and rather go through regions of high values of $f$. 
Hence, the main idea (already proposed and studied in previous works) is to define a Riemannian metric on $\mathbb R^d$ that is conformal to the Euclidean metric and penalizes regions of low values of $f$. 

Here, $f$ is unknown, and only random points are available, so $f$ and its derivative must be estimated.

The paper provides no theory - it stitches together the two ideas described above: finding shortest paths relative to $f$ and estimating $f$ and its derivative from the data. However, the ideas used in the algorithms and the numerical experiments are very much appreciated by the reviewers, who all (except for one, who does raise important theoretical questions) are in favor of accepting the paper.

Hence, I also advocate for accepting the paper. Please include the additional material and discussions in the final version of the paper.
Decision: Accept (poster)