Learning Distances from Data with Normalizing Flows and Score Matching

Rebuttal to Reviewer nn54

We thank the reviewer for their careful and detailed reading of our work, and for the helpful questions and clarifications.

---

Novelty of the Algorithms

The algorithms used by the authors are already present in the literature.

While 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.

Algorithm 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.

---

Metric Tensor Notation

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$?

It 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$.

---

Scope of the Analysis and Manifold Assumptions

More generally, I am struggling to understand in which space the analysis in the paper is conducted.

We 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$.

---

Clarification on “Distance Between Points in a Distribution”

Does ‘define distance between points in a distribution’ mean ‘define distance between points sampled from a distribution’?

Yes, that is correct. Thank you — we will update the phrasing accordingly to avoid ambiguity.

---

When Are Geodesics Distance-Minimizing?

Do you know in which Riemannian manifold the geodesics are always distance minimizers?

This 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).

---

We 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.

Rebuttal to Reviewer XLew

Thank you for your positive assessment of our work and for your thoughtful suggestions.

---

Preliminary Results on Real-World Dataset

Please 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.

---

Computational Efficiency Compared to Graph-Based Methods

Please 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.

Graph-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.

---

Assumptions Underlying Theoretical Results

The theorems provided in the paper, particularly the derivation of the geodesic equations and their reparameterization, rely on standard smoothness assumptions:

• The density function $p(x)$ is assumed to be differentiable and non-zero in the region of interest.
• The score function $s(x) = \nabla \log p(x)$ is assumed to be Lipschitz continuous.
• The conformal metric defined by $p(x)$ is smooth.

These 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.

---

Thank you again for your encouraging review and support.

Rebuttal to Reviewer AqTp

We thank the reviewer for their thoughtful and constructive feedback. Below, we address each of the main points and questions raised.

---

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?

The tractability of the relaxation method primarily depends on the complexity of the data, rather than the ambient dimension. For example:

• Uniform distributions: Our method correctly converges to straight lines even with very coarse discretization.
• 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.

More 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.

---

How would the method perform in small and simple image datasets such as MNIST or CIFAR?

We 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.

As 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.

As 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:

We 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.

---

I assume the geodesic equation for Fermat metrics is a known result (reference?), and that only the reparameterization is new.

Yes, 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.

---

I suggest adding a plot showing an example of geodesics recovered on a dataset by the different methods.

Agreed—we will add such a figure to the main text to complement Figure 4, comparing several methods on the same dataset.

---

A plot with approximated geodesics with various discretization levels would be interesting.

We will add such a figure to the appendix to illustrate convergence behavior and practical trade-offs.

---

Minor Comments

• Figures 1 and 4: Will ensure they are clearly referenced.
• Repeated sentence on p3: Will be removed.
• Fig. 10: Will standardize y-axis ranges for better visual comparison.

---

We 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.

Rebuttal to Reviewer EAyN

We thank the reviewer for their detailed and constructive comments.

---

Simplicity of Experiments

Please see our response to reviewer AqTp for results on an experiment on MNIST.

---

Missing Recent Work

Thank you for pointing out our oversight. We will add both references to the discussion of related work.

---

Time Complexity and Scaling

Scaling 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.

For 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.

The 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.

---

Fermat Distance Scaling in Non-Gaussian Distributions

You 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$.

By 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)).

A 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.

Thank you for encouraging us to clarify and strengthen this argument — we will update the text and appendix accordingly.

---

Speed and Initialization of Algorithm 1

Please 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).

---

Computation of Ground Truth Geodesic Distance (Eq. 13)

We 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.

---

Absence of Graphs for NF-Based Relaxation

You 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.

---

Hybrid Training of NFs with Score Matching

We 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.

---

Line 434: Unifying Normalizing Flows and Score Models

By 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.

---

Line 102: Numerically Stable Relaxation Method

We 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.

---

We thank the reviewer again for their thoughtful feedback and for helping us sharpen the clarity and scope of the work.