Estimating Riemannian Metric with Noise-Contaminated Intrinsic Distance
Authors: Jiaming Qiu
Xiongtao Dai
Venue: NeurIPS 2023 poster (2023)
Keywords: metric learning
manifold learning
local metric
dissimilarity
geometry

Abstract: We extend metric learning by studying the Riemannian manifold structure of the underlying data space induced by similarity measures between data points. The key quantity of interest here is the Riemannian metric, which characterizes the Riemannian geometry and defines straight lines and derivatives on the manifold. Being able to estimate the Riemannian metric allows us to gain insights into the underlying manifold and compute geometric features such as the geodesic curves. We model the observed similarity measures as noisy responses generated from a function of the intrinsic geodesic distance between data points. A new local regression approach is proposed to learn the Riemannian metric tensor and its derivatives based on a Taylor expansion for the squared geodesic distances, accommodating different types of data such as continuous, binary, or comparative responses. We develop theoretical foundation for our method by deriving the rates of convergence for the asymptotic bias and variance of the estimated metric tensor. The proposed method is shown to be versatile in simulation studies and real data applications involving taxi trip time in New York City and MNIST digits.

Estimating Riemannian Metric with Noise-Contaminater Intrinsic Distance
Jiaming Qiu
Xiongtao Dai
Fred Hutchinson Cancer Center Friuu@fredheee or.
Division of Biostatistics, School of Public Health University of California, Berkeley xda i @berreleeee edu
Abstract
We extend metric learning by studying the Riemannian manifold structure of the underlying data space induced by similarity measures between data points. The key quantity of interest here is the Riemannian metric, which characterizes the Riemannian geometry and defines straight lines and derivatives on the manifold. Being able to estimate the Riemannian metric allows us to gain insights into the underlying manifold and compute geometric features such as the geodesic curves. We model the observed similarity measures as noisy responses generated from a function of the intrinsic geodesic distance between data points. A new local regression approach is proposed to learn the Riemannian metric tensor and its derivatives based on a Taylor expansion for the squared geodesic distances, accommodating different types of data such continuous, binary, or comparative responses. We develop theoretical foundation for our method by deriving the rates of convergence for the asymptotic bias and variance of the estimated metric tensor. The proposed method is shown to be versatile in simulation studies and real data applications involving taxi trip time in New York City and MNIST digits.
Introduction
With increasing data complexity, it has received increasing attention for capturing the geometric structure of the data space characterized by Riemannian manifold, which consists of two vital parts: the coordinate system and the Riemannian metric. Manifold learning attempts to find low-dimensional representations as coordinates for the high-dimensional data [23, 27, 7]. Yet existing methods generally assume the local geometry is given, e.g., by the Fuclidean distance between ambient data points, leading to non-isometric methods that distort the data geometry in the embedded space. the other hand, many methods have been proposed to consider the Riemannian metric, including parametric estimation [17, 16, 15]; multi-metric learning [12]; pushforward via the Laplacian operator [21]; pull-back via Jacobian of generative models [1, 2]; Bayesian approach based are length of Gaussian process [14]; and dynamics of evolving probability when time-indexed samples are available [24]. It is also possible to recover the manifold in an abstract setting with only similarity measures among data points [9].
We propose to estimate the Riemannian metric utilizing similarity measures between data points when data coordinates are available. Suppose that data are generated from an unknown Riemannian manifold, while the Euclidean distance between the coordinates may not reflect the underlying geometry. We model the observed similarity measures between data points as noise-contaminated intrinsic distances, which are then used to characterize the intrinsic geometry on the Riemannian manifold. The targeted Riemannian metric is estimated in a data-driven fashion, which enables path finding via geodesics and calculus on the manifold in which the data reside.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

This problem is closely related to metric learning [3, 26] where a Mahalanobis distance is commonly used to obtain the best distance for classification [11] and clustering [29] tasks. While global metric is often the focus of earlier works, multiple local metrics [10, 28, 22, 5] are found to be useful because they better capture the data space geometry. Despite the resemblance, our target is to learn the Riemannian metric instead of the distance metric, which fundamentally differentiates our proposal from metric learning. We emphasize the data space geometry rather than obtaining a metric optimal for subsequent supervised learning tasks. The additional smoothness of the Riemannian metric tensor field allows analysis of finer structures, while the coefficient matrix for the Mahalanobis distance is locally constant for distance metric learning.
To formulate the problem, let (M, G) be a Riemannian manifold with Riemannian metric and induced geodesic distance dist (-.-) which measures the true intrinsic difference between points. Knowing the coordinates of data points zo, 21 E M, we identify each point via its coordinates (a tuple of real numbers). The noisy measures of the intrinsic distance between data points are referred to as similarity measures (equivalently dissimilarity). The response modeled flexibly, and we consider the following common scenarios: (i) noisy distance, where y = disl (Iy, #1)² + e for error where y 0 if the two points 2o, 21 considered similar and y 1 otherwise, and (iii) relative comparison. where a triplet of points (æo, 21, z2) are given and y = 1 if 2n is more similar to 2 than to 22 and y 0 otherwise. The binary similarity measurement is common in computer vision [e.g. 6]. while the relative comparison could be useful for perceptional tasks and recommendation system [e.g. 25, 4]. We aim to estimate the Riemannian metric G and its derivatives using the coordinates and similarity measures among the data points.
(ii) similarity/dissimilant
The major contribution of this paper is threefold. First, we formulate a framework for probabilistic modeling of similarity measures among data on manifold via intrinsic distances. Such a framework also justifies the distance metric learning from geometrical perspective and unveils its approximation error for the first time as we know. Second, a theoretical foundation is developed for the proposed method including asymptotic consistency. Last and most importantly, the proposed method provides a geometric interpretation for the structure of the data space induced by the similarity measures, as demonstrated in the numerical examples that include taxi travel and an MNIST digit application.
2 Background in Riemannian Geometry
For brevity, metric now refers to Riemannian metric while distance metric is always spelled out. Throughout the paper, M denotes a d-dimensional manifold endowed with a coordinate chart (U, ()), where o U Rd maps a point p E U C .M the manifold to its coordinate (p(p) (g' (p)..... .., 40t (p) E Rd. Without loss of generality, we identify point by its coor- dinate as (p'. ...... p²). suppressing yp for the coordinate chart. Upper-script Roman letters denote the components of coordinate, e.g., p² is the i-th entry in the coordinate of the point p, and 7 is the i-th component function of a curve y : R 3 Ja, b] -> .M when expressed on chart U. The tangent space T,MM is vector space consisting of velocities of the form v = 7(()) where 7 is any curve satisfying (0) = p. The coordinate chart induces basis the tangent space T TpM,, as 0
0/0z"), Ip for i = I....,d, d, so that a tangent vector e E Tp.M is represented as v = Σi v'O, for some o² E IR suppressing the subscript p in the basis. We adopt the Einstein summation convention unless otherwise specified, namely e'O, denotes Σ v'O,. where common pairs of upper- and lower-indices denotes a summation from to d [see e.g., 18, 1991
The Riemannian metric G on a d-dimensional manifold JM is smooth tensor field acting on the tangent vectors. At any p M, G(p): TpMx x T.M - RR is symmetric bi-linear tensor/function satisfying G(p)(v, v) 2 0 for any v E TpM and G()) = 0 if and only if v = O. On a chart sp, the metric is represented as a d-by-d positive definite matrix that quantifies the distance traveled along infinitesimal changes in the coordinates. With an abuse of notation, the chart rep- resentation of G is given by the matrix-valuce function p + G(p)
(Gg( E Radd for P E M, so the distance traveled by ry at time t for a duration of dt is (G)(
The intrinsic distance induced by G, or the geodesic distance, is computed as dist (p. inf, To VEisassa Giy(a(t))á' (f)de (t)dt for two points p. q on the manifold M, where infimum (eyisssaGG(aa())))
is taken over any curve a [0, 1] + M connecting p q.

---

A geodesic curve (or simply geodesic) is a smooth curve y R. D [a, b] + M satisfying the geodesic equations, represented on a coordinate chart as
(1) + (00)))))) o >(()) = 0, for i, 1,...,d,
(2.1)
where over-dots represent derivative w.r.t. t; 1 G (a.G 0 Ga 0;Gū;) are the Christoffel symbols at p; and Ghi is the (ks l)-emmenn of G Solving (2.1) with initial conditions' produces geodesic that traces out the generalization of a straight line on the manifold, preserving travel direction with acceleration, and also locally the shortest path.
Considering the shortest path ry connecting p to and applying Taylor's expansion at t = 0, we obtain dist (p, x)² x Eisiyea Giz (g) - p') - p²), showing the connection between the geodesic distance the Mahalanobis distance. Our estimation method is inspired by this approximation, and we will discuss the higher-order terms shortly which unveil finer structure of the manifold.
3 Local Regression for Similarity Measurements
3.1 Probabilistic Modeling for Similarity Measurements
Suppose that we observe N independent triplets (Yu. Xuo, Xu1), ), I.....N. N. Here, the Xuj are locations on the manifold identified with their coordinates uj:.... (X,
E R = 1,2, and Y, are noisy similarity measures of the proximity of (Xwo, Xx1) in terms of the intrinsic geodesic distance dist (+:+) on M. To account for different structures of the similarity measures, we model the response in a fashion analogous to generalized linear models. For Xwo, Xu1 lying in small neighborhood Up M of target location p E M, the similarity measure Y, is modeled as
E [Yu|Xxwo X.1) -ー (disl (Xao, Xxn)²),
(3.1)
where g is given link function that relates the conditional expectation the squared distance. Example 3.1. We describe below three common scenarios modeled by the general framework (3.1).
1. Continuous response being the squared geodesic distance contaminated with noise:
Yu = dist (Xuo X.1)² (p)
(3.2)
where E...... .....,,f. Eu are i.idd mean zero random variables, and M R is smooth positive function determining the magnitude of noise near the target point p. This model will be applied to model trip time as noisy measure of cost to travel between locations.
2. Binary (dis)similarity response:
P (Y, = 1 Nwo Xw1) = logit (dist (X Xw1)² (p))
(3.3)
for some smooth function & : : R. where logit(p) = log (m/ (1 - H)). / (0, 1) is the logit function. This models the case when there are latent labels for Xay and Y. is a measure of whether their labels are in common not. The function ft (p) in (3.3) describes the homogeneity of the latent labels for points in a small neighborhood of p. The latent labels could have intrinsic variation even if measures
made for the same data points x = X = Xa1, and the strength of which is captured by h. (p).
3. Binary relative comparison response, where we extend our model for triplets of points (Xao, Xu1, Xw2), where Y, stands for whether X,00 more similar to Xu1 than to Xw2:
IP (Y, = 1|x X.1, Xw2) = logit- -1 (dist (X X,z)² dist (Nav, Xua)²), (3.4)
that the relative comparison Y, reflects the comparison of squared distances.
See Section $3 of the Supplement for details about solving in practice.

---

3.2 Local Approximation of Squared Distances
It is the squared distance that determines the responses in our model (3.1), which justified by the following local approximation. Proposition 3.1 provides a Taylor expansion for the squared geodesic distance between two geodesics with same starting point but different initial velocities (see Figure 1 for visualization). It is the key tool to estimate model (3.1) through local regression. Furthermore, (3.5) characterize the error of approximating geodesic distance with Mahalanobis distance. For a point p on the Riemannian manifold M, let exp, Tp.. - M denote the exponential map defined by (tv) = rt) where Y is a geodesic starting from at time 0 with initial velocity 7(())
M. For notational simplicity, we suppress the dependency on p in geometrie quantities (e.g., the metric G' is understood to be evaluated at p). For i == ,,....d, d. denote 8' = (1) = (t) (0) as the difference in coordinate after travel of time allng 7.
Proposition 3.1 (spread of geodesics, coordinated), Let p E M and v. w E TpM be two tangent vectors at p. On local coordinate chart, the squared geodesic distance between two geodesics no(t) (tu) and (t) (tp, satisfies, as t - 0,
dist (t)( (t))² 16-10 Gig + + (05%) - o58%) 12 + 0(t*), where for j, k. m ,,....d.
(3.5)
% (() p', of = 71 (1) - p', and _/ - 6 - 81. i.e., %. oj are differences in i-th coordinates of (t) and (t) to the origin respectively, and %/ = % - 8 is the coordinate difference between (t) and (t);
. Gij and Th are the elements of the metric and Christoffel symbols at p. respectively.
Proof. See Section S5 in the Supplement.
To the RHS of (3.5), the first term is the quadratic term in distance metric learning. The second term the result of coordinate representation of geodesics. vanishes under the normal coordinate where the Christoffel symbols are zero. It inspires the use of local regression to estimate the metric tensor and the Christoffel symbols. For Xuo, Xu1 in a small neighborhood of p, write the linear predictor as
7a == g(0) + 5th -
8 0 88 + 62,0-1 (Dta@%
5th 6/ 3022
(3.6)
igh: where 8n = X 5" = -p', and %0,0-1 1 Sig 5,., for j, k, 1, 1....., d, and ,,...,.. N.. The intercept term (00) is included for capturing intrinsic variation (e.g., hp)) in (3.3)) and can otherwise be dropped. The link function connects the linear predictor to the conditional mean via thu == -- (nw) s F. (Y, Xwo, Xui) as indicated by (3.1) and (3.5), where nu is seen function of the coefficients ((0), / and a(2)
a function of the intercept ((0) and coefficintss (6 [(2)
Therefore, upon the specification of a loss function Q R R - (0) U UR++ and non-negative weights wi...., WN the minimizers
201)
arg min gto g 3 a. a
Σ Q (Ya: (h) Was
(3.7)
u=1
subject to
3 = 3
(1)
jike for jor, 1,, L....,d
1,
(3.8)
used to estimate the metric tensor and Christoffel symbols, obtaining
3()
N = 1g
8
(3.9)
ij
where G is the matrix inverse of & satisfying Gd Gd1 = (()=))- The symmetry constraints (3.8) the result of the symmetries in the metric tensor and Christoffel symbols, and are enforced by optimizing over only the lower triangular indices < < j S without constraints. Asymptotic
Seee e.g., [19] pp. 131-133 for normal coordinate. See [20] and Proposition S5.1 of the Supplement for boordinateinnvariant invariant version of Proposition 3.1.

---

results provide the cositive-definitenn of the metric estimate, as will be shown in Proposition 4.1. To weigh the pairs of endpoints according to their proximity to the target location p. we apply kernel weights specified by
d II K
X (
Kin--") v²)
W/ h-2d
K
(3.10)
h
h
x=1
for some h > O and non-negative kernel function K. The bandwidth controls the bias-variance trade-off of the estimated Riemannian metric tensor and its derivatives.
Altering the link function and the loss function Q in (3.7) enables flexible local regression estimation for models in Example 3.1.
Example 3.2. Consider the following loss functions for estimating the metrie tensors and the Christof- fel symbols when data are drawn from model (3.2)-(3.4), respectively.
1. Continuous noisy response: use squared loss Q (y, (a) = (y - 14] with being the identity link function so pau ?/-
2. Binary (dis)similarity response: use log-likelihood of Bernoulli random variable
Q (y. (1) = y log μ + (1 - y) log (1 -p),
(3.11)
and the logit link, so nuu logit (m)) The model becomes a local logistic regression.
3. Binary relative comparison response: apply the same loss function (3.11) and logit link in the previous scenario, but here we formulate the linear predictor based on dist (Xuo: Xu2)² dist (Xuo.X.1) X,1)² 1711 - 1u22 and
Nu = g (Na nu2) Locally, the difference in squared distances is approximated by
(3.12)
(8) 1 -5, 2,,0-2000000 /1g + (8t (8an - 8% 81))
Mal 7u2
(5t,0 2 (5in0h - 512222)) 1(22)
(3.13)
for 81,2 = X22 -p' and 8 8'22 8ap- 1,.... d. Here 7a1 and 7/22 are constructed in analogy (3.6) using (Xuo, Xu1) and (X Xu2) pair respectively.
Examples in Section 5 will further illustrate the proposed method in those scenarios. Besides the models listed, other choices for the link g and loss function Q can also be considered under this local regression framework [8], accommodating a wide variety of data. To efficiently estimate the metric on the entire manifold M, we apply procedure based on discretization and post-smoothing. as detailed in Section S3 of the Supplement. In short, kernel smoothing (weighted average) of estimated
component functions G.j over some grid points provides smooth tensor field, which ease the burden of re-estimation for every points on the manifold.
Bias and Variance of the Estimated Metric Tensor
This subsection provides asymptotic justification for model (3.2) with (Ya|/Xuo X.1)
dist XX1) under the squared loss Q(p, y) (μ y) and the identity link g(m) μ.. The estimator we analyzed here fits local quadratic regression without intercept and approximates the squared distance by a simplified form of (3.6):
dist (Xao, Xu1)² se nu == di Hu.0-1hhp-10 (11)
(4.1)
for u = 1....,N. N. Given a suitable order of the indices i.j for vectorization, we rewrite the formulation into a matrix form. Denote the local design matrix and regression coefficients
D. (6 1
8 8,0-18 10-18_0-1
8 80,0-10 (2,01000011)
B (0(1) ..... 3( ..... 801) T
5

---

so that the linear predictor 7 = =. Further, write D = (Di...... +Di.....Dw) Y (Y,,.... YN) and W diag (Wi, .... WN). with weights Wa, specified in (3.10). The objective function in (3.7) becomes - D W -W--DD
- 3D and the minimizer is ß (D"WD) "DWYY, for which we will analyze the bias and variance.
To characterize the asymptotic bias and variance of the estimator, we assume the following conditions satisfied in neighborhood of the target p. These conditions are standard and analogous those assumed in local regression setting [8].
(A1) The joint density of endpoints (Xmo, X.1) positive and continuously differentiable.
(A2) The functions Gúj: D are C²-smooth for i, k = ,,...,d.
(A3) The kernel K in weights (3.10) is symmetric, continuous, and has bounded support.
(A4) sup, var (Yu Xu1) < oc.
Proposition 4.1. Under (A4), bias (β|x) Op (h²) var (Qx) = Op ( as h - 0 and N42+2d oc, where X {(Xu,, X.1) Xu1) is the collection of observed endpoints.
Proof. See Section S6 in the Supplement.
The local approximation (4.1) similar local polynomial estimation of the second derivative of a 2d-variate squared geodesic distance function, explaining the order of f in the rates for bias and variance.
The positive definite constraints (as in distance metric learning) are dropped based on the asymptotic properties, thus the estimated metric could fail to be positive definite at some locations. This is alleviated by sufficient data and post-smoothing that averaging neighboring estimation.
Simulation
illustrate the proposed method using simulated data with different types of responses as described in Example 3.1. We study whether the proposed method well estimates Riemannian geometric quantities, including the metric tensor, geodesics, and Christoffel symbols. Additional details includ- ing assessment stability and positive-definiteness are included in Section S4 of the Supplement. Subsection S4.2 also provides an ellipsoid example for manifold of hon-constan curvature.
5.1 Unit Sphere
The usual arc-length/great circle distance on the d-dimensional unit sphere is induced by the round metric, which is expressed under the stereographic projection coordinate (x*...., xx)) G.j
(1+ E1 ((i=)) for i,j 1..... d. Under the additive model (3.2) in Example 3.1,
we considered either noiseless or noisy responses by setting a (pp) 0 or (p) > 0 respectively.
Experiments were preformed with d 2 and the finding are summarized in Figure 5.1. For continuous responses, the left panel of Figure la visualizes the true and estimated metric tensors via cost ellipses (S2.1) and the right panel shows the corresponding geodesics by solving the geodesic equations (2.1) with true and estimated Christoffel symbols. The metrics and the geodesics were well estimated under the continuous response model with or without additive errors, where the estimates overlap with the
truth. Figure 5.1c evaluates the relative estimation errors G G | and [ [f-r]]
|[]]
w.r.t. the Frobenius norm (S2.2) for data from the continuous model (3.2).
For binary responses under model (3.3), Figure 5.1b visualizes the data where the background color illustrates A. Figure 5. and left panel of Figure 5. suggest that the intercept and metric were reasonably estimated, while the geodesics are slightly away from the truth (Figure 1. right). This indicates that the binary model has higher complexity and less information provided by the binary response (see also Figure S4. lb in the Supplement).

---

metric tensor
geodesics
0.6
0.3
+
0.0
-2
-0.3
-2
0
-0.6
x1
-2
0.50
0.75
1.00
1.25
1.50
x1
intercept
10 2030
truth
estimated with additive error
method
estimated with error
estimated with binary response
binary response
(a) Estimated and true metric tensors using ellipses representation (left) and the geodesic curves (right) starting from (1, 0) with unit initial veloe- ities pointing 1-12 o' cloock directions.
(b) Binary responses, with segments showing comparisons. background colored to ft
metric
Christoffel Symbol
estimated intercept
-1
-2
2
1
2
-2
-1
0
x1
x1
relative error
relative
relative error
5200
050'0
$60'
000
error
05'0
52'D
100
520 0500
(c) Relative errors in term of Frobenius norm of the estimated tensors for (d) Errors estimating with bi-
the continuous response model (3.2) with additive
nary responses (3.3).
Figure 5.1: Simulation results for 2-dimensional sphere under stereographic projection coordinate.
5.2 Relative Comparison on the Double Spirals
set of 7 x 10 points on R² were generated around two spirals, corresponding to two latent classes A and B (e.g., green points in Figure 5.2a are from latent class A). We compare neighboring points (Xwo, Xx1, Xx,) to generate relative comparison response Y, as follows. For u = 1,..., N, Y,, 1 if X no, X.1 belong the same latent class and Xwo, Xwa belong to different classes; otherwise Y 0. Figure 5.2b shows portion of the = 6, 965, 312 comparison generated, where the hollow circles in the middle of each wedge correspond to X.o-
Here, contrast the two latent classes induces the intrinsic distance, so the distance is larger across the supports of the two classes and smaller within single support. Therefore, the resulting metric tensor should reflect less cost while moving along the tangential direction of the spirals compared to perpendicular directions. Estimates
drawn under model (3.4) by minimizing the objective (3.7) with the link function (3.12) and the local approximation (3.13).
The estimated metric shown in Figure 5.2c is consistent with the interpretation of the intrinsic distance and metric induced by the class membership discussed above. Meanwhile, the estimated geodesic curve unveils the hidden circular structure of the data support shown in Figure 5.2d.
7

---

response
latent class
(d) Geodesics starting from (1.0) with initial velocity pointing to clock under estimated metric.
class
(a) The points and their la- (b) Relative comparison. tent classes
Cireles show Xuo, branch
(c) Estimated metric
to X.1, Xa2.
Figure 5.2: Relative comparison on double spirals. Gray curves (solid and dashed) in the background represent the approximate support of the two latent classes.
6 New York City Taxi Trip Duration
We study the geometry induced by taxi travel time in New York City (NYC) during weekday morning rush hours. New York City Taxi and Limousine Commission (TLC) Trip Record Data was accessed on May 1st, 2022 to obtained business day morning taxi trip records including GPS coordinates for pickup/dropoff locations as (X.o, X.a) and trip duration as Y,. Estimation to the travel time metric was drawn under model (3.2) with Q n)) = (y 12)² and g(u) = μ..
Figure 6. shows the estimated metric for taxi travel time. The background color shows the Frobenius norm of the metric tensor, where larger values mean that longer travel time is required to pass through that location. Trips through midtown Manhattan and the financial district were estimated be the most costly during rush hours, which is coherent to the fact that these are the city's primary business districts. Moreover, the cost ellipses represent the cost in time to travel unit distance along different directions. This suggests that in Manhattan, it takes longer to drive along the east-west direction (narrower streets) compared to the north-south direction (wider avenues).
Geodesic curves in Figure 6. b show where a 15-minutes taxi ride leads to starting from the Empire State Building. Each geodesic curve corresponds to of 12 starting directions (1-12 o' ccooc). Note that we apply a continuous Riemannian manifold approximation to the city, so the geodesic curves provide approximations
the shortest paths between locations and need conform the road network. Travel appears to be faster in lower Manhattan than in midtown Manhattan. The spread of the geodesics differs along different directions, indicating the existence of non-constant curvature on the manifold and advocating for estimating the Riemannian metric tensor field instead of applying a single global distance metric.
7 High-dimensional Data: An Example with MNIST
The curse of dimensionality is challenging for nonparametric methods applied on data sources like images and audios. However, it is often found that apparent high-dimensional data actually lie close to some low-dimensional manifold, which is utilized by manifold learning literature to produce reasonable low-dimensional coordinate representations. The proposed method can be applied the resulting low-dimensional coordinates to recover the high-dimensional geometry even implicitly defined, demonstrated in the following MNIST example.
We embed images MNIST to a 2-dimensional space via tSNE [13]. Similarity between the objects was computed by the sum of the Wasserstein distance between images" and the indicator of whether the underlying digits are different (1) or not (0). The goal is to infer the geometry of the embedded
e hhrps://wwwl nnyc htwp://heehchaaattttl---r--roooe page. Data format changed after our download.
After rescaling, see Subsection S4.5 of the Supplement.

---

40.85'N
Frobenius
time (minutos)
40.8
b4
6075
0.3
74 00 W3.05 73.00'W
(b) Geodesics correspond 15-minute taxi rides from the Empire State Building heading to 1-12 of c'oloc.
(a) Estimated metric tensors for trip duration: cost ellipses and Frobenius norm (background color).
Figure 6.1: New York taxi travel time during morning rush hours.
data induced by this similarity measures. Estimation drawn under model (3.2) with squared loss Q (4) = (y - 2 and the identity link. An intereept term (8(0) in (3.6) is used to capture intrinsic variation partially due to the non-injective dimensional reduction.
The geodesics estimated from our method tend to minimize the number of switches between labels. For example, the geodesic
Figure 7.1 remains "4" (Ist row of panel (b)) throughout, while the straight line on the (SNE chart translates to path of images switching between "4" and "9" (2nd row panel (b)); similar phenomenon occurs for geodesics B and C (corresponding 7th and 12th rows in (b)). Also, our estimated geodesics produce reasonable transition and reside in the space of digits, while unrestricted optimal transport (3rd, 8th, and 13th rows of panel (b)) could produce unrecognizable intermediate images. Our estimated geodesic is faithful to the geometric interpretation that geodesic is locally the shortest path. Moreover, our proposal provide sensible image transition. especially compared to the shortest paths along neighhorhood graph in the embedded space (dotted blue lines in (a) and 4, 9, rows in (b)) or kNN graphs in the original Fuclidean space 28 (5, 10, 15 rows in (b)). This could be useful for complex data interpolation that preserves geometry.
8 Discussion
We present a novel framework for inferring the data geometry based on pairwise similarity measures. Our framework targets data lying low-dimensional manifold where observations expected to be dense near the locations where we wish to estimate the metric. However, these assumptions general enough for our method to be applied to manifold data with high ambient dimension in combination with manifold embedding tools. Context-specific interpretation of the geometrical notions, e.g., for Riemannian metric and geodesics, has been demonstrated in the taxi travel and MNIST digit examples. Our method could inspire additional compelling applications to other topics such cognition and perception research, where psychometric similarity measures are commonly made. Moreover, further development including (but not limited to) handling sparse data and ensuring positivity in the Riemannian metric estimation is needed, discussed in Section SI of the Supplement.

---

9
9
4
4
4
Y
3
8
8
3
3
3
9 9
9 9
3
3
3
3
3
8
S S 3
3 3 3
3
3
S
C
3
3
7
S
5
8
S
3
3 3
3
8
2
S S
5 5
5 5
5 5
22
S S 5
8
8
8 9
f
5 S
?
(b) Image transitions corresponding to A, B. and C in
(a). Every 5 rows correspond to set of paths sharing
-2
2
the same pair of starting and ending images, where the
ISNE.1
1st - 5th row correspond the estimated geodesics, straight lines on the chart, the optimal transport (path not shown in (a)), the paths along neighborhood graph on the chart, and along kNN graph in the Euclidean space (path not shown (a)). respectively.
(a) The estimated geodesic curves (solid red), shortest path along graph (dotted blue), and straight lines on the chart (dashed black).
Figure 7.1: Geometry induced by a sum of Wasserstein distance and same-digit-or-not indicator.
10

---

Acknowledgments and Disclosure of Funding
The authors would like thank Dr. Zhengyuan Zhu and the reviewers for the helpful feedback and discussion. The research was partially supported by NSF grant DMS-2329879.
References
[1] Georgios Arvanitidis, Lars Kai Hansen, and Søren Hauberg. Latent space oddity: On the curvature of deep generative models. In International Conference on Learning Representations, February 2018.
[2] Georgios Arvanitidis, Søren Hauberg, and Bernhard Schölkopf. Geometrically enriched latent spaces, August 2020.
[3] Aurélien Bellet, Amaury Habrard, and Mare Sebban. Metric learning. Synthesis Lectures on Artificial Intelligence and Machine Learning, 9()::-151,, January 2015.
[4] Adam Berenzweig, Beth Logan, Daniel P.W. Ellis, and Brian Whitman. A large-scale evaluation of acoustic and subjective music-similarity measures. Computer Music Journal, 28(2):63-76. June 2004.
[5] Shuo Chen, Lei Luo. Jian Yang, Chen Gong, Jun Li, and Heng Huang. Curvilinear distance metric learning. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[6] S. Chopra, R. Hadsell, and Y. LeCun. Learning similarity metric discriminatively, with application to face verification. In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR 05), volume I, pages 539-546 vol. 1, June 2005.
[7] Ronald Coifman and Stéphane Lafon. Diffusion maps. Applied and Computational Harmonic Analysis, 21(1):5-30, July 2006.
[8] Jianqing Fan and I. Gijbels. Local Polynomial Modelling and Its Applications. Number 66 in Monographs on Statistics and Applied Probability. Chapman & Hall, London New York, Ist ed edition, 1996.
[9] Charles Fefferman, Sergei Ivanov, Matti Lassas, and Hariharan Narayanan. Reconstruction of a Riemannian manifold from noisy intrinsic distances. SIAM Journal Mathematics of Data Science, (3):770-808, September 2020.
[10] Andrea Frome, Yoram Singer, Fei Sha, and Jitendra Malik. Learning globally-consistent local distance functions for shape-based image retrieval and classification. In 2007 IEEE 11th International Conference on Computer Vision, pages 1-8, October 2007.
[11] T. Hastie and R. Tibshirani. Discriminant adaptive nearest neighbor classification. JEEE Transactions Pattern Analysis and Machine Intelligence, 18(6):607-616, June 1996.
[12] Søren Hauberg, Oren Freifeld, and Michael J. Black. A geometric take metric learning. In F. Pereira, C. J. C. Burges, L.. Bottou, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 25, pages 2034-2032. Curran Associates, Inc., 2012.
[13] Geoffrey E Hinton and Sam Roweis. Stochastic Neighbor Embedding. In Advances in Neural Information Processing Systems, volume 15. MIT Press, 2002.
[14] Martin Jørgensen and Soren Hauberg. Isometric Gaussian process latent variable model for dissimilarity data. In Proceedings of the 38th International Conference Machine Learning, pages 5127-5136. PMLR, July 2021.
[15] Tam Le and Marco Cuturi. Unsupervised Riemannian metric learning for histograms using Aitchison transformations. In Proceedings of the 32nd International Conference Machine Learning, pages 2002-2011. PMLR, June 2015.
[16] G. Lebanon. Metric learning for text documents. IEFE Transactions on Pottern Analysis and Machine Intelligence, 2844):497-508 April 2006.

---

[17] Guy Lebanon. Learning Riemannian metrics. In Proceedings of the Nineteenth Conference on Uncertainty in Artificial Intelligence, UAI'03, pages 362-369, San Francisco, CA, USA, 2002. Morgan Kaufmann Publishers Inc.
[18] John M. Lee. Introduction Smooth Manifolds. Number 218 in Graduate Texts in Mathematics. Springer, New York; London, 2nd ed edition, 2013.
[19] John M. Lee. Introduction to Riemannian Manifolds. Springer Berlin Heidelberg, New York, NY, 2018.
[20] Wolfgang Meyer. Toponogov's theorem and applications. Lecture Notes, Trieste, 1989.
[21] Dominique Perraul-Joncar and Marina Meila. Non-linear dimensionality reduction: Riemannian metric estimation and the problem of geometric discovery. May 2013.
[22] Deva Ramanan and Simon Baker. Local distance functions: A taxonomy, new algorithms, and an evaluation. IEEE Transactions on Putlern Analysis and Machine Intelligence, 33(4):794-806, April 2011.
[23] Sam T. Roweis and Lawrence K. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science, 290(5500):2323 2326, December 2000.
[24] Christopher Scarvelis and Justin Solomon. Riemannian Metric Learning via Optimal Transport. In The Eleventh International Conference on Learning Representations, February 2023.
[25] Matthew Schultz and Thorsten Joachims. Learning distance metric from relative comparisons. In Advances in Neural Information Processing Systems, volume 16. MIT Press, 2003.
[26] Juan Luis Suárez, Salvador García, and Francisco Herrera. A tutorial on distance metric learn- ing: Mathematical foundations, algorithms, experimental analysis, prospects and challenges. Neurocomputing, 425:300-311 -20000022,222 February 2021.
[27] Joshua B. Tenenbaum, Vin de Silva, and John C. Langford. A global geometric framework for nonlinear dimensionality reduction. Science, 290(5500) December 2000.
[28] Kilian Q Weinberger and Lawrence K Saul. Distance metric learning for large margin nearest neighbor classification. Journal of machine learning research, 10(2), 2009.
[29] Eric Xing, Michael Jordan, Stuart Russell, and Andrew Ng. Distance metric learning with application to clustering with iide-nnformtiio. In S. Becker, S. Thrun, and K. Obermayer, editors, Advances in Neural Information Processing Systems, volume 15. MIT Press, 2002.
12


Reviews:
[{"review_id": "FUDHYNOdq4", "reviewer": "Reviewer_E2so", "summary": "This paper proposes a method to estimate the Riemannian metric of data space when coordinate representations of each data point and some noisy similarity measurements among data points are provided. The similarity measurements of different types, such as noise-contaminated distances, similarity/dissimilarity labels, and comparative similarity labels, are probabilistically modeled as functions of intrinsic distances between data points. The squared geodesic distances are approximated to be linear to the Riemannian metric and Christoffel's symbol via Taylor's expansion, leading the metric estimation problem to be like a maximum likelihood estimation using a generalized linear model with each entry of the Riemannian metric and the Christoffel's symbol as its parameters. Asymptotic convergence rates of the bias and variance of the estimator are derived for the case where noisy distances are given as the similarity measure. Experiments are performed using some simulated data, New York taxi trip duration data, and MNIST data sets to demonstrate the proposed estimator's benefits in capturing the underlying geometry of the data space.\n\n**Post rebuttal** I have increased my score from 4 to 5. Lingering concerns involve providing more practical applications and comparisons to other methods.", "questions": "* When extracting the Riemannian metric from binary similarity labels or binary relative comparison labels, there may not be any ground truth Riemannian metrics. Is it then appropriate to call the proposed framework an 'estimation'?\n* How can the proposed probabilistic models be justified? Distance scale may significantly affect the results according to the current choice of models in (3.3) and (3.4). Also, different options of probabilistic models would result in other Riemannian metrics.\n* Is there any guarantee that the local squared distance in (3.5) is always positive when the estimation is performed based on (3.7)? This may also be related to the point that the estimated metric is not guaranteed positive-definite. I wonder if solving for the case (3.11) and (3.12) can make the estimated squared distances, i.e., $\\eta_u$s, negative.\n* There are no explanations for obtaining estimates for Christoffel's symbols at arbitrary points, i.e., $\\hat{\\Gamma}\\circ \\gamma(t)$, which are required for obtaining geodesics.\n* Regarding the post-smoothing of the obtained Riemannian metrics explained in the supplement, the weighted averaging of Riemannian metrics seems strange. For example, we cannot simply add tangent vectors on different tangent spaces but should apply parallel transport for them to be on the same tangent space. \n\n**[Minor comments]**\n* It must be explicitly stated in the main text the location of the supplement containing the proof of each proposition.\n* Since the proposed approximations and estimation methods are valid only locally, it would be better to explain how the methods can be applied globally in the main text, not in the supplement.\n", "limitations": "The authors have adequately addressed the limitations of the paper. For other possible limitations, please refer to the weaknesses and questions above.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "* The paper is in general well-written.\n* The paper provides some interesting ideas to characterize the Riemannian geometry of data space by utilizing noisy similarity measurements among data points, such as continuous noisy distances, binary similarity labels, and binary relative comparison labels, in a unified framework.\n* The proposed method to estimate the Riemannian metric is quite simple, seems original, and shows reasonable experimental results.\n* This paper provides complete proof of a proposition for the local approximation of the geodesic distance and that for the asymptotic convergence analysis of the estimator.\n", "weaknesses": "* The considered problem setting seems to be quite restrictive. It would be rare that a problem provides both low-dimensional coordinate representations and meaningful intrinsic distances among data points. Obtaining each of them has been an important research topic for decades.\n* The experiments need to provide practical applications of the proposed method. For example, obtaining geodesics in the experiments is usually done to demonstrate the validity of the estimated Riemannian metrics but not for further use. Suggesting more practical uses for real data based on the estimated Riemannian metric, computed geodesics, or other subsequent geometric quantities such as lengths or volumes, would significantly benefit the community.\n* This paper lacks any comparison with other Riemannian metric estimation methods. A direct comparison would be possible to the method in Perrault-Joncas et al., which considers a mapping $f$ from the data submanifold embedded in the ambient space (endowed with an ambient space metric) to a coordinate chart and estimates the push-forward metric in the coordinate chart via the mapping $f$.\n* Experiments in high-dimensional settings are lacking. How the method behaves according to the dimension would be valuable information. A worry here is that local regions to estimate the Riemannian metrics would increase exponentially with respect to the dimension.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nFor the weakness pointed out regarding restrictive setup, application, comparison with other methods, and dimensionality, we kindly ask the reviewer to refer to our general response points 1--3. See our general response point 5 for compatibility of the metric.\n\nThe reviewer is correct that the proposed method (and most if not all local polynomial based approaches) do not scale well with increased dimension. That is why this work is based on a low-dim manifold assumption as discussed in general point 1. As we demonstrated in the MNIST example, it successfully preserve the geometry of high-dim image data with low-dim spaces.\n\n## How can the proposed probabilistic models be justified?\nMetric estimated under a constant multiple of the distances (scaling) will only differ up to multiplication of a constant. This is a special case of conformality. For example, one can follow the definition of the Christoffel symbol (line 88) to see that it is invariant to different distance scales. In other words, the different distance scales would act like unit conversion and not affect the geometry features like angles, shapes, and the geodesic curves (but its length would differ).\nEmploying different probabilistic models (specifically $g$ and $Q$) may lead to different estimated metrics, but the estimated metric will change only if the derivative of the link at 0 changes by the chain rule. Like that in the generalized linear models, this is a necessary modeling component which would be determined by data type and domain knowledge. Our contribution focuses on the possibility of estimating the metric under such framework with a spread of geodesics, for which simple models as in (3.2) – (3.4) were considered.\n\n## Positivity\nWhile the estimated Riemannian metric in a finite sample is not guaranteed to be positive definite, our asymptotic theory (proposition 4.1) shows that the estimated metric will be positive definite with probability tending to 1 given sufficient observations (near the point of estimation). In our experiments we rarely had estimates that were not positive definite and we did not encounter serious non-positive issues as shown in the last part of section S4.5 (page 15) and figure S4.10 of the Supplement. One could in principle adopt constraint optimization to enforce positive definiteness as that in the metric learning, but for simplicity we chose to employ an unrestricted algorithm that is more computationally efficient. \nModel (3.2) is general in that it can include both positive and negative distances, but in practice only positive distances are observed and our model is applied on these more intuitive setups. One could switch the link function and/or the loss function in a similar manner as in local quasi-likelihood [(Fan, Heckman, and Wand 1995)](https://www.tandfonline.com/doi/abs/10.1080/01621459.1995.10476496) to ensure strictly non-negative distance. But again, we chose to present the simplest model here.\n\n## There are no explanations for obtaining estimates for Christoffel's symbols at arbitrary points.\nFollowing (3.6) – (3.9), the Christoffel symbols can be estimated simultaneously with the metric, and analogous post-smoothing steps (see section S3) also apply. Post-smoothing is adopted mainly to speed up the geodesic computation. So that we do not need to re-estimate the tensors at every point as requested by the ODE solver, which is time consuming. We also notice the ODE solver (and the resulting geodesic curves) benefited from the added smoothness especially for the cases with binary similarity measures.\n\n##  The weighted averaging of Riemannian metrics seems strange in the post-smoothing step.\nSince we are working on a chart, post-smoothing is applied on the component functions instead of on the (coordinate-free) tensor directly, so comparison of tangent vectors on different tangent spaces is not involved. More specifically, under the local coordinate chart $(x^1, \\dots, x^d)$, the estimated metric tensor is fully spelled out as $G_{ij} dx^i dx^j$. We post-smooth the estimate of the (continuous) component functions $G_{ij}, i, j = 1, \\dots, d$. Implicitly, we assume the coordinate chart contains the data domain, so that transition maps (to bridge different local coordinate charts) are not involved.\n\n## Minor comments\n- *It must be explicitly stated in the main text the location of the supplement containing the proof of each proposition.*  \n\nWe will explicitly state this in the camera-ready version if accepted. Thank you for pointing this out.\n\n- *Since the proposed approximations and estimation methods are valid only locally, it would be better to explain how the methods can be applied globally in the main text, not in the supplement.*\n\nWe will add a short sketch of the post-smoothing and refer to the Supplement in the main text."}, {"review_id": "bzwvmIansp", "reviewer": "Reviewer_5uut", "summary": "This paper develops a theory for estimating the Riemannian metric tensor for a given set of observations and some addition information. This additional information includes a (noisy) measure of similarity between the given points in a pairwise fashion. Examples of this information includes the geodesic distance, or a binary response about the similarity of types, or a binary response about relative comparison. \n\nThe formulation relies on a formula for the intrinsic distance between corresponding points on two geodesics shot from the same base p.  This formula approximates the said distance using the Euclidean chords, Riemannian metric tensor, and the Christoffel symbols. The first term computes the Riemannian metric between the shooting vectors and the second term accounts for the curvature. The paper seizes on this linear approximation and sets up a regression problem for estimating these matrices from the given data. In this way, it estimates the metric tensors and the Christoffel symbol at point p. The paper further develops the estimation theory (bias and variance) for the metric tensor estimation. \n\nThis theory is demonstrated using several experiments, some based on simulated data and some on real data sets. The simulated data helps validate the method for a known geometry (sphere, double spiral). The real data experiments involve learning metrics for the taxi travel times and MNIST images embedded in R^2 using tSNE. \n", "questions": "\n\n-- I believe that the equation connecting geodesic distance to the Mahalanobis distance is the same kind that is used to derive Cramer-Rao bound on estimation error using the Fisher-Rao Riemannian metric. It would be interesting to make that connection if there is one. \n\n\n-- What are the practical situations where the pairwise geodesic distances between points on a manifold are given, and one does not know the Riemannian metric? I understand the taxi example but for the MNIST example one has to use an embedding which is somewhat of an arbitrary choice. ", "limitations": "\nThe authors have not explicitly address the limitations. ", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "-- The problem of learning Riemannian metric from the data is quite important and challenging. \n\n-- The paper formulates this problem in a novel and interesting way and provides analytical insights, rather than the current systems that apply deep neural networks to such problems. \n\n-- The paper goes on to develop estimation theory for the metric tensor. There is some interesting statistical contributions here. ", "weaknesses": "\n-- The simulation experiments seem to involve manifolds with constant curvature (unless the curvature changes along the spiral). Perhaps the authors can try their method for ellipsoids or some variable curvature manifolds. \n\n-- How does the method work if the sampled points on the manifold are sparse? I feel that the linear regression model derived here requires the manifold points to be close otherwise the errors will start piling up. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "For the practical problems with pairwise distance but no Riemannian metric, we kindly refer to our general response point 2. See also our general point 4 for limitations.\n\n## The simulation experiments seem to only involve manifolds with constant curvature.\nOur theoretical foundation allows for general non-constant curvature and this is the backbone that makes the method work for practical applications as demonstrated. The proposed method applies to manifolds with variable curvature like an ellipsoid, and we can showcase this via an ellipsoid example if accepted and time permitted.\n\n## How does the method work if the sampled points on the manifold are sparse?\nSufficiently dense data is indeed necessary for the proposal, which could be demanding for high-dimensional data due to the curse of dimensionality. However, data often exhibits the manifold phenomenon, namely data intrinsically lie close to a low-dimensional manifold (see also our general point 1). Thus as illustrated in the MNIST example, we can apply the proposed method after a dimension reduction step such as tSNE. The resulting representations tend to be dense since they lie in a low-dimensional space. This can substantially alleviate the curse of dimensionality, and the dense local neighborhood requirement will more likely hold true.\n\n## Connection to Fisher-Rao metric\nWhen working with a statistical manifold, namely a parametric family where each point is a distribution, the Fisher-Rao metric is the canonical Riemannian metric which enjoys nice statistical interpretation. Its component matrix is the Fisher information, the inverse of which becomes the Cramer–Rao lower bound. Amari has shown that the Kullback–Leibler divergence (like a distance) is a function whose second differential produces the Riemannian metric (e.g., [Murray Rice 1993](https://books.google.com/books?id=ZBa7F9LrDrMC&pg=PA76#v=onepage&q&f=false)). Though, we are uncertain whether there would be a direct connection to our work estimating the metric tensors, while the Fisher-Rao's are based on parametric family."}, {"review_id": "b38HchORge", "reviewer": "Reviewer_9sCn", "summary": "The paper aims to extend metric and manifold learning by learning Riemannian metrics from functions of the observed data that are not related to an embedding space metric as is the usual case. Assuming there is an underlying Riemannian metric structure and that the observed dissimilarity is a known function of this structure, the authors show how a Riemannian metric can be learned from dense data. The method is tested on trip time data in New York and MNIST represented on a 2d manifold. \n", "questions": "Convincing counterarguments to the weaknesses listed above.\n\n- do you require the charts to be normal? If so, that restricts the set of metrics you can possibly learn, right? I.e., it could be that a true underlying Riemannian structure was not orthonormal with respect to the chosen chart", "limitations": "the limitations of the method could be described more precisely", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- methods for estimating a Riemannian metric from dissimilarity measures are derived\n- in principle, it can be a good idea relying on other distances than embedding space distances to learn geometric structure", "weaknesses": "- learning Riemannian metrics has a long history in the literature, although not in this exact form\n- I am unsure of the usefulness of the method:\n-- are there any guarantees that a Riemannian metric that is compatible with the observed structure exists? I believe if the objective is actually a (sufficiently smooth) distance there are existence results, but it would be nice with arguments for the more general cases in sec 3.1\n-- it is assumed that a chart is known. One can always find a mapping to a lower-dimensional subspace such as is done in the MNIST example, but whether the data is actually 2d is unknown, and the uncertainty in the chart estimation is not taken into account\n-- in practice, the result would likely be a Riemannian metric approximation in a chart that is estimated with high uncertainty and with strong assumptions on the dimensionality. The data would likely not be dense in higher dimensions, and so the estimated structure would very likely be a poor approximation. This can then be used for e.g. geodesic interpolation, but whether this is actually useful is not clear to me\n-- what happens if the model is misspecified? The link functions are assumed to be known a priori, right?\n- the developed estimation techniques have merit, but the fact that a Riemannian metric can be recovered from the metric (distances) is not a new result\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We acknowledge the existence of related literature learning Riemannian metrics in lines 22-26 and will include more as suggested by reviewers. We hope the proposed framework can shad new light on this topic. See also our general response point 3. See our general point 4 for limitations. \n\n## usefulness of the method\nAs we discussed in the general point 5, existence results for such compatible metric do exist. [Fefferman et al. 2020](https://doi.org/10.1137/19M126829X) also discussed a manifold reconstruction problem similar to our additive error case (eq. (3.2)). Their discussion focuses on abstract metric spaces with no pre-specified coordinates, while our setup is simpler as we assume that the coordinate chart is readily available.\n\nIndeed lower-dimensional subspaces are always obtainable while the question is how faithful they would be to present the original data. We consider this a modeling choice while specifying the coordinate chart as a step in, e.g., tSNE. In fact, assessing the uncertainty of the proposed models (3.1) could also potentially provide insight to the quality of dimension reduction: for example, a large uncertainty may suggest poor coordinate representations.\n\nThe link function is pre-specified and the model can be misspecified. Though, our focus is the connection in (3.1) and (3.5), therefore only those simpler models (3.2) – (3.4) (as commonly seen in generalized linear models) are included. As we discussed in the general point 1, more flexible models are possible under the proposed framework, and that the effect of the link function is expected to be small.\n\n## Do you need normal coordinates?\nNo, the chart need not be normal. See line 140 of the main text. The Christoffel symbols are coordinate dependent and vanish under normal coordinates. Our method does not assume this and the Christoffel symbols are estimated."}, {"review_id": "9fjxlB49SV", "reviewer": "Reviewer_uGdp", "summary": "The paper presents a novel mechanism for learning a Riemannian metric from distance observations. This is important for applications where relative observations are available (e.g. \"objects x1 and x2 are different, while x2 and x3 are similar\"), such as perception studies. The approach is based on local regression from which a Riemannian metric is deduced from a Taylor expansion of the geodesic distance. Limited experimental results demonstrate feasibility.", "questions": "* What does the $\\circ$ operator denote in Eq. 2.1? As I read the equation, all quantities are scalars, so I do not think that this denotes a Hadamard product, but this is nonetheless my best guess.\n* I do not quite understand Eq. 3.2: Are there no further restrictions in $\\epsilon_n$? If so does that not mean that $P(Y < 0) > 0$, which is rather odd for distance observations which must be strictly non-negative.\n* There has been quite extensive work on pull-back metrics in autoencoder-like models (see e.g. \"Latent Space Oddity: on the Curvature of Deep Generative Models\", Arvanitidis et al., ICLR 2018). These can be viewed as smooth interpolations of local linear regressions since the metric is formed by Jacobian matrices. Do you see links between this line of work on the present paper?", "limitations": "The main paper only makes a brief mention of the limitations, while a more extensive discussion is given in the supplements. I appreciate the supplementary discussion but would encourage the authors to move this to the main paper to increase openness regarding limitations.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "* The paper approaches an important question that is highly understudied in the machine-learning community.\n* The proposed Riemannian metric estimator is new and sufficiently simple to be practically useful.\n* It is neat that the estimator also comes with a direct estimate of the Cristoffel symbols as these can otherwise be tedious to compute.\n* The examples given throughout the paper are instructive and provide a nice assistance to the reader.\n* The illustrative examples are instructive.", "weaknesses": "* The approach seems to be difficult to scale to higher dimensional data. I don't think this is a significant problem as the current application areas where distance observations are available often revolve around low-dimensional observations, e.g. in perception studies in psychology.\n* The computational procedure seems to rely on a discretization of the input space, which will only work in low-dimensional cases.\n* The kernel smoothing used to produce the Riemannian metric (eq. 3.10) is rather ad hoc. It would have been nice if this was more closely tied to the local regression. That being said, the approach basically follows established procedures, see e.g. [10].\n\n### Related work\n* I am not fond of the phrasing of lines 37-39, which suggests that the proposed work is the first time that Riemannian metrics are learned. The cited work of Lebanon [14] and Hauberg et al. [10] follows the same philosophy as the present paper, so I think a rephrasing would be appropriate.\n* The authors are not the first to learn metrics from observations of distances. I think it would be good to at least acknowledge the vast literature on multi-dimensional scaling, which does not learn a Riemannian metric but learns from distances. Perhaps closer to the present work is the paper \"Isometric Gaussian Process Latent Variable Model for Dissimilarity Data\" (Jørgensen et al., ICML 2021) which learns a low-dimensional manifold and its Riemannian metric from noisy distance observations. The taken approach is quite dissimilar from the present paper, but the ambition is similar.\n\n### Minor things\n* Typo in line 219: \"without or without\" --- I suppose this should be \"with or without\".", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "For concerns regarding scalability to high-dim data, we kindly refer the reviewer to our general point 1. We also thank the reviewer for pointing out many related works, as discussed in our general point 3, we will add citations accordingly. We also discuss the positivity question in the general point 8. We also fix the typo in line 219 (for a camera-ready version hopefully), thank you for the careful inspection.\n\nThe $\\circ$ operator in eq. (2.1) is function composition, i.e., $f \\circ g (t) = f(g(t))$. Adopted to avoid too many brackets especially in our case with many sub/superscripts.\n\n## Ad hoc kernel smoothing in eq. (3.10)\nThe objective (3.7) and weights in equation (3.10) are motivated by the Taylor expansion (3.5) similarly to that in the local polynomial regression (e.g., section 5.4.3 of \\[6\\]). The weights (3.10) are relatively standard, treating the $(\\delta_{u0}^1, \\dots, \\delta_{u0}^d, \\delta_{u1}^1, \\dots, \\delta_{u1}^d)$ as a $2d$-array. This method has a solid theoretical foundation as established in our Proposition 4.1. One key difference between our proposal to that in \\[10\\] is that the weights are applied differently. We estimate the Riemannian metric tensor directly and the weights (3.10) are involved in the objective function (3.7), while the method in \\[10\\] constructs the metric tensor by smoothing multiple local distance metric (see eq. (6) in \\[10\\]), in a similar manner as our post-smoothing step (section S3 in the Supplement).\n\n## Potentially misleading claim in lines 37-39\nLines 37-39 (and this paragraph spanning lines 34-42) emphasize that it is the *Riemannian metric* that we are targeting which is different to the *distance metric* as in the classic metric learning. It does not imply that this is the first time Riemannian metrics are learned, as we have referred to multiple references (as in lines 22-26) including those you kindly pointed out.  \n\n## Positivity\nWhile the estimated Riemannian metric in a finite sample is not guaranteed to be positive definite, our asymptotic theory (proposition 4.1) shows that the estimated metric will be positive definite with probability tending to 1 given sufficient observations (near the point of estimation). In our experiments we rarely had estimates that were not positive definite and we did not encounter serious non-positive issues as shown in the last part of section S4.5 (page 15) and figure S4.10 of the Supplement. One could in principle adopt constraint optimization to enforce positive definiteness as that in the metric learning, but for simplicity we chose to employ an unrestricted algorithm that is more computationally efficient. \nModel (3.2) is general in that it can include both positive and negative distances, but in practice only positive distances are observed and our model is applied on these more intuitive setups. One could switch the link function and/or the loss function in a similar manner as in local quasi-likelihood [(Fan, Heckman, and Wand 1995)](https://www.tandfonline.com/doi/abs/10.1080/01621459.1995.10476496) to ensure strictly non-negative distance. But again, we chose to present the simplest model here.\n\n## Link to works on pull-back metrics in autoencoder-like models\nThe pull-back approach tackles a similar problem recovering data space geometry. It utilizes the Jacobian matrix of some smooth mapping $f: \\mathcal Z \\to \\mathcal X$ to pull the metric on $\\mathcal X$ back to $\\mathcal Z$. Typically the mapping $f$ is a generator/encoder, while $\\mathcal Z$ and $\\mathcal X$ are some latent space and the data space respectively. It is usually assumed that the metric $\\mathcal M_x$ of input space is readily available, e.g., assumed to be Euclidean ([Arvanitidis, Hansen, and Hauberg 2018](https://openreview.net/forum?id=SJzRZ-WCZ)), or estimated via Riemannian metric learning ([Arvanitidis, Hauberg, and Schölkopf 2020](http://arxiv.org/abs/2008.00565)). \nOur proposal can provide the metric of latent space directly based on similarity measures in the data space in a manner similar to our MNIST example. The major difference is that our method treats the low-dim embedding as a coordinate chart, which is a stricter assumption: the pull-back metric only requires smooth mappings between the data and latent space (e.g., generator or encoder)."}]

Author rebuttal:
We truly appreciate reviewers' careful read of our manuscript and helpful comments. We summarized common weakness/questions and provide below our general response. See our replies to individual reviews for specific responses. 

## 1. Restrictive problem/model setting
*The reviewers question that our settings could be strict: a. low-dimensionality and meaningful distance, b. probabilistic formulation in eq (3.2) – (3.4), and c. a pre-specified link function.*  

1. Our method and calculation are designed for data drawn from a manifold with low intrinsic dimensions. The observed similarity measures are generated based on intrinsic geodesic distance, as in eq (3.1). The low-dimensional manifold assumption ([Bengio et al. 2013](https://doi.org/10.1109/TPAMI.2013.50)) is commonly satisfied by real-world image and audio datasets because of the manifold phenomenon, and also satisfied in perception studies in psychology as pointed out by a reviewer. The manifold assumption enables us to apply nonlinear dimension reduction and avoid the curse of dimensionality, even though the raw data, like images, live in a high-dimensional ambient space. 
2. The probability model is proposed as a principled guide for method development and also for theory. The specific models enlisted in Example 3.1 and Example 3.2 are given as commonly encountered scenarios that can be handled by our framework; the proposed method is not limited to these models. 
3. The choices of link function $g$ and loss function $Q$ are flexible and can accommodate a variety of data generating mechanisms similar to what generalized linear model (GLM) can handle. We expect the effect of the link function to be small because the estimation is done locally and only the derivative of the link function at the origin matters.

## 2. Lack of practical application/importance
We demonstrated two applications from transportation and computer vision. We expect our method to also be widely applicable in perception studies but the authors are not familiar with open data there.  

Problems do exist where pairwise distances are easily obtainable but the driving geometry is not. Other than our taxi and MNIST examples, correlations among fMRI signals (track recurrent coactivation of neurons) are used to analyze brain functional connectivity ([van den Heuvel and Hulshoff Pol 2010](https://doi.org/10.1016/j.euroneuro.2010.03.008)). Another (remotely) related topic is travel time tomography, where the internal structure of media (e.g., organ, earth) is estimated by the travel time of waves (e.g., ultrasound, seismic). They also utilize pairwise measures to capture underlying structure. 

In many cases, the data space is non-Euclidean, thus a Riemannian metric is needed to capture the intrinsic geometry. Yet, finding an appropriate one can be a challenging task but of great interest, as it could lead to more accurate similarity measures, better clustering algorithms, and improved recognition systems. We refer to literature in our main text and the general point 3 and references therein.

## 3. Existing methods learning Riemannian metric
*There is a large body of existing literature learning Riemannian metric, and there is no comparison to existing methods such as that in \[18\].*

We will include additional citations to acknowledge more recent literature utilizing latent variable models, including those based on pull-back ([Arvanitidis, Hauberg, and Schölkopf 2020](http://arxiv.org/abs/2008.00565)) and/or GP-LVM ([Jørgensen and Hauberg 2021](https://proceedings.mlr.press/v139/jorgensen21a.html)). 
Meanwhile, stemming from distance metric learning, our work focuses on connecting the similarity measures to the Riemannian metric and its higher order derivatives (as in our eq. (3.1) and (3.5)). Admittedly there are plenty of rooms for further development including benchmarking to comparable methods, which cannot be exhausted in a single paper.

Perrault-Joncas et al. (2013)’s method ( \[18\] ) targets a more limited scenario in that it starts with pairwise extrinsic Euclidean distances without noise, while our method handles general distance metric with noise. Their method can thus not be applied to most of the numerical demonstrations in our work except for the MNIST digit. There is no code implementation of their work available which makes it difficult for comparisons.

## 4. Lack of open discussion on limitations
We do recognize several limitations and discuss those in section S1. We appreciate the reviewers’ interest and a deep dive into our Supplement. The limitation section is put here because of the page limit for the main text. We will add a more apparent remark in the main text linking to them.

## 5. What are we estimating?
*Can we always find a Riemannian metric compatible with the similarity measure, so that we are actually “estimating” some well-defined quantities, especially for binary measures?*

It is not always guaranteed that a (distance) metric space is a Riemannian manifold, let along when the measurements are binary. Analogously, in order for an underlying Riemannian metric to exist, the similarity measure needs to be induced by the geodesic distance on some Riemannian manifold. Thus our method is based on the manifold assumption (Bengio et al. 2013) and that the observed similarity measures are induced by some Riemannian metric, for which we are estimating. We are not aware of a certain answer for the binary situation. However, an abstract metric space can be approximated by a Riemannian manifold (theorem 1 of \[7\]). A later paper ([Fefferman, et al. 2020](https://doi.org/10.1137/19M126829X)) also discussed a manifold reconstruction problem similar to our additive error case (eq. (3.2)).

## 6. Positivity
*The estimated metric tensors might fail to be positive-definite and model (3.2) allows negative distance.*

Due to length limit of a single comment, see the Positivity section of our response to reviewers E2so and uGdp (identical).

Meta-review: This paper tackles the estimation of Riemannian metrics on latent spaces from noisy observations of distance data, either through the distances themselves or weaker observations such as relative distance differences. This problem is important and challenging, and the submitted paper presents a novel and (in this context) simple way to estimate the underlying Riemannian metric. This is valuable!

While two reviewers are enthusiastic about the paper, the other two are borderline, and during rebuttal, the authors have provided a number of clarifications, as well as some promises of further experiments. Their final version should make sure to include the clarifications and promised additional experiments, as well as an explicit and open discussion of the paper's limitations.
Decision: Accept (poster)