Estimating Riemannian Metric with Noise-Contaminated Intrinsic Distance
Authors: Jiaming Qiu
Xiongtao Dai
Venue: Submitted to ICLR 2023 (2023)

Abstract: We extend metric learning by studying the Riemannian manifold structure of the underlying data space induced by dissimilarity 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 dissimilarity 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. Our framework is general and accommodates different types of responses, whether they are continuous, binary, or comparative, extending the existing works which consider a single type of response at a time. 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 a real data application involving taxi trip time in New York City.

Under review as a conference paper at ICLR 2023
ESTIMATING RIEMANNIAN METRIC WITH NOISE-CONTAMINATED INTRINSIC DISTANCE
Anonymous authors Paper under double-blind review
ABSTRACT
We extend metric learning by studying the Riemannian manifold structure of the underlying data space induced by dissimilarity 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 to gain insights into the underlying manifold and compute geometric features such as the geodesic curves. We model the observed dissimilarity 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. Our framework is general and accommodates different types of responses, whether they are continuous, binary, or comparative, extending the existing works which consider a single type of response at a time. 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
The estimation of distance metric, also known as metric learning, has attracted great interest since its introduction for classification (Hastie & Tibshirani, 1996) and clustering (Xing et al., 2002). A global Mahalanobis distance is commonly used to obtain the best distance for discriminating two classes (Xing et al., 2002; Weinberger & Saul, 2009). While a global metric is often the focus of carlier works, multiple local metrics (Frome et al., 2007: Weinberger & Saul, 2009; Ramanan & Baker, 2011; Chen et al., 2019) are found to be useful because they better capture the data space geometry. There is a great body of work on distance metric learning; see, e.g., Bellet et al. (2015); Suárez et al. (2021) for recent reviews.
Metric learning is intimately connected with learning on Riemannian manifolds. Hauberg et al. (2012) connects multi-metric learning to learning the geometric structure of Riemannian manifold, and advocates its benefits in regression and dimensional reduction tasks. Lebanon (2002; 2006); Le & Cuturi (2015) discuss Riemannian metric learning by utilizing a parametric family of metric, and demonstrate applications in text and image classification. Like in these works, our target is to learn the Riemannian metric instead of the distance metric, which fundamentally differentiates our approach from most existing works in metric learning. We focus on the nonparametric estimation of the data geometry as quantified by the Riemannian metric tensor. Contrary to distance metric learning, where the coefficient matrix for the Mahalanobis distance is constant in a neighborhood, the Riemannian metric a smooth tensor field that allows analysis of finer structures. Our emphasis is in inference, namely learning how differences in the response measure are explained by specific differences in the predictor coordinates, rather than obtaining a metric optimal for a supervised learning task.
A related field is manifold learning which attempts to find low-dimensional nonlinear representa- tions of apparent high-dimensional data sampled from underlying manifold (Roweis & Saul, 2000: Tenenbaum et al., 2000; Coifman & Lafon, 2006). Those embedding methods generally start by assuming the local geometry is given, e.g. by the Fuclidean distance between ambient data

---

Under review as a conference paper at ICLR 2023
points. Not all existing methods are isometric, so the geometry obtained this way can be distorted. Perraul-Joncas & Meila (2013) uses the Laplacian operator to obtain pushforward metric for the low- dimensional representations. Instead of specifying the geometry of the ambient space, our focus is to learn the geometry from noisy measures of intrinsic distances. Fefferman et al. (2020) discusses an abstract setting of this task, while our work proposes a practical estimation of the Riemannian metric tensor when coordinates are also available, and we show that our approach is numerically sound.
We suppose that data are generated from an unknown Riemannian manifold, and we have available the coordinates of the data objects. The Euclidean distance between the coordinates may not reflect the underlying geometry. Instead, we assume that we further observe similarity measures between objects, modeled as noise-contaminated intrinsic distances, that are used to characterize the intrinsic geometry on the Riemannian manifold. The targeted Riemannian metric is estimated in a data- driven fashion, which enables estimating geodesics (straight lines and locally shortest paths) and performing calculus on the manifold.
To formulate the problem, let GM, be a Riemannian manifold with Riemannian metric G, and dist (+-) be the geodesic distance induced by G which measures the true intrinsic difference be- tween points. The coordinates of data points To, 21 E M are assumed known, identifying each point via tuple of real numbers. Also observed are noisy measurements of the intrinsic distance between data points, which we refer to as similarity measurements (equivalently dissimilarity). The response is modeled flexibly, and we consider the following common scenarios: (i) noisy distance, where y dist (wo, x1)²² + e for error e, (ii) similarity/dissiiiarrtt where y 0 if the two points To. #1 are considered similar and
otherwise, and (iii) relative comparison, where triplet of points (Ig, 21. 22) are given and y if ro is more similar 21 than to z2 and y otherwise. The binary similarity measurement is common in computer vision (e.g. Chopra et al., 2005), while the relative comparison could be useful for perceptional tasks and recommendation system (e.g- Schultz& & Joachims, 2003; Berenzweig al., 2004). We aim to estimate the Riemannian metric G and its derivatives using the coordinates and similarity measures among the data points.
The major contribution of this paper is threefold. First, we formulate a framework for probabilis- tic modeling of similarity measurements among data on manifold via intrinsic distances. Based on Taylor expansion for the spread of geodesie curves in differential geometry, the local regression procedure successfully estimates the Riemannian metric and its derivatives. 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 measurements, as demonstrated in the numerical examples that in- clude a taxi travel and an MNIST digit application.
BACKGROUND IN RIEMANNIAN GEOMETRY
For brevity, metric now refers to Riemannian metric while distance metric is always spelled out. Throughout the paper, MM denotes a d-dimensional manifold endowed with coordinate chart (U. p). where 3
(p(p) (4 (p),.... .,. ddd + Rd maps a point p E U C M on the manifold to its coordinate (p)) E Rd.. Without loss of generality,
identify point by its coordi- nate as (p' .......)) )²) suppressing yp for the coordinate chart. Upper-script Roman letters denote the components of a coordinate, e-g. p² is the i-th entry in the coordinate of the point p. and y is the i-th component function of a curve IR プ [a, b] - .M when expressed on chart U. The tangent space Tr.M a vector space consisting of velocities of the form v = 7((0) where 7 is any curve satisfying (0) = p. The coordinate chart induces a basis on the tangent space T,MM, as 0 = 0/0z"), I, for i = 1.....,, d, so that a tangent vector v E TpM is represented as v = Σi= aa for some e² E R, suppressing the subscript p in the basis. We adopt the Einstein summation convention unless other- wise specified, namely u v'O, denotes Σd= v'O,, where common pairs of upper- and lower-indices denotes summation from to d (see e.g., Lee, 2013, 1p.18-19)
The Riemannian metric G on a d-dimensional manifold M is a smooth tensor field acting on the tangent vectors. At any p M. G(p): : : T.. xT.M
TpMM + R. is symmetric bi-linear tensor/function satisfying (())(v, v) 2 0 for any 0 E T and ((p) v) = 0 and only if v = O. On a chart sp, the metric is represented as d-by-d positive definite matrix that quantifies the distance traveled
2

---

Under review as a conference paper at ICLR 2023
along infinitesimal changes in the coordinates. With an abuse of notation, the chart representation of is given by the matrix-valued function p + G(p) = (r)(f))) E Rdxd for p E M, so the distance traveled by at time & for duration of dt is [G (1))))))))))))))))) The intrinsic distance induced by G, or the geodesic distance, is computed as
dist (p.9) inf S/¹
(yylatt)d"t))d)))ttt ISaisd
Σ
(2.1)
for two points p. 9 on the manifold M, where infimum is taken over any curve [0.1] - M connecting p to
A geodesic curve (or simply geodesic) is smooth curve Y R J [a + M satisfying the geodesic equations, represented on a coordinate chart as
(t) + (())))(((()) 7(t) 0, for i, I,...,d,
(2.2)
where over-dots represent derivative w.r.t. L; I = 2 G (a,G,) + - are the Christof- fel symbols p. and G is the (k, -elemmen of G Solving (2.2) with initial conditions' produces geodesic that traces out the generalization of a straight line on the manifold, preserving travel direc- tion with acceleration, and is also locally the shortest path.
Considering the shortest path y connecting p to and applying Taylor's expansion at t = 0, we obtain dist (p:4)* q)² R Eisijsa Ge (g) - ('r - p²), showing the connection between the geodesic distance and a quadratic form analogous to the Mahalanobis distance. Our estimation method is based on 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, Xwo, Xu1), 1, ,, N. Here, the Xuj are locations on the manifold identified with their coordinates (Naj...., X R j 1,2, and Y, are noisy similarity measures of the proximity of (Xuo, X.1) in terms of the intrinsic geodesic distance dist (-.)) on M. To account for different structures of the similarity measurements, we model the response in a fashion analogous generalized linear models. For Xuo, X.1 lying in small neighborhood Up C M of a target location p E M, the similarity measure Y,, is modeled as
E (YulXuo: X..) = g- (dist (Xao, X,.)²).
(3.1)
where g is given link function that relates the conditional expectation to the squared distance. Éxample 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:
Y. dist (Xwo, X.1)² (p))n,
(3.2)
where E .......6 en are i.i.d. mean zero random variables, and a M - R+ is a smooth positive function determining the magnitude of noise near the target point p. This model will be applied to model trip time noisy measure of cost to travel between locations.
2. Binary sdis)similarity response:
P PY = 1|x Xw1) = logit (dist (Xuo, Xu1)² (p))
(3.3)
for some smooth function : .M- + R, where logit(n) = log (H/ (1 (11-))), H)), / E (0. 1) is
the logit function. This models the case when there are latent labels for Xus (c.g., digit v.s. 9) and Y, is a measure of whether their labels are in common or not. The function (p) in (3.3) describes the homogeneity of the latent labels for points in small neighborhood of The latent labels could have intrinsic variation even if measurements are made for the same data points x X..o X.1, and the strength of which is captured by h. (p).
See Appendix B for details about solving it in practice.
3

---

Under review as a conference paper at ICLR 2023
3. Binary relative comparison response, where extend model for triplets of points (Xuo: Xul, Xu2), where Yu stands for whether Xuo is more similar to Xu1 than to Xu2:
TP (Yn = 1|x Xu1, Xx2) = logit." (disll (Xuo, Xu2)² 2 dist (Xuo, Xaa)), (3.4)
that the relative comparison Y, reflects the comparison of squared distances.
3.2 LOCAL APPROXIMATION OF SQUARED DISTANCES
We develop local approximation for the squared distance the key tool estimate our model (3.1) through local regression. Proposition 3.1 provides a Taylor expansion for the squared geodesic dis- tance between two geodesics with same starting point but different initial velocities (see Figure D.1 for visualization). For a point p on the Riemannian manifold M, let expp TpMM - .M denote the exponential map defined by (tv) = -() where yy is a geodesic starting from at time 0 with initial velocity 7(()) v E p.M For notational simplicity, we suppress the dependency on p in geometric quantities (c.g.. the metric G is understood to be evaluated at p). For i I....,dd dd, denote 8 = (t) = (t) (0) as the difference in coordinate after a travel of time along 7 Proposition 3.1 (spread of geodesics, coordinated). Let p E M and v, TpM be two tangent vectors at p. On a local coordinate chart, the squared geodesic distance between two geodesics (n(() (xp)(lv) and 71 (() expp(Ixc) satisfies, as t
dist (ro(t), (1) 60-18-1 + 6- (058) 078)) where for k,j,,,,, ......d.
G2 O(t*)
(3.5)
%8 ()(() p', of = =1) p², and 6/ ob - 8i, i.e., 86- of differences in i-th coordinates of (t()) and 71 (t) to the origin p. respectively, and 88_ 53 5j is the coordinate difference between 70 (4) and (f);
Gij and T2 are the elements of the metric and Christoffel symbols at p. respectively.
To the RHS of (3.5), the first term is the quadratic term in distance metric learning. The second term is the result of coordinate representation of geodesics. It 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 Xwo, Xa11 in a small neighborhood of p, write the linear predictor as
Tu = (00) + 60 L 1 + 00,0-1 ( -8 62) 3
8
(3.6)
a function of the intercept ((0) and coefficients 8 ((2)))) where 82o Xio - pp', 8 = = p'
((2)
and 6 6 8% - 5, for i, k, 1.... d, and
I....,.. N. The intercept term g(0) is included for capturing (p) in (3.3) under that scenario and can otherwise be dropped from the model. The link function connects the linear predictor to the conditional mean via μ4, := g (mw) x
E (Y.| X..) as indicated by (3.1) and (3.5), where Ma is seen as a function of the coefficients
((00,, aQ)) and (2(2) Therefore, upon the specification of a loss function Q R R - (0) R and non-negative weights w1..... WN, the minimizers
N Σ u=1
(9
30)) ((2) Pij
arg min
Q (Yu, Hu) War
(3.7)
subject to
8
3(2) jik
for i. k. I,....d,
(3.8)
are used to estimate the metric tensor and Christoffel symbols, obtaining
G = 3()
:j = -,22
(3.9)
SSe e.g., Lee (2018) pp. 131-133 for normal coordinate. See Meyer (1989) and Proposition D.1 for coordinate-invariant version of Proposition 3.1.
4

---

Under review as a conference paper at ICLR 2023
where Ghi is the matrix inverse of G satisfying G Gk) (z=1)) The symmetry constraints (3.8) are the result of the symmetries in the metric tensor and Christoffel symbols, and are enforced by optimizing over only the lower triangular indices 1 S <i S d without constraints. Asymptotic sults guarantees the positive-definitness 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
W/ h 2d
II K
d
p') (K h
X
K
ul
PD) h
(3.10)
for some h > O and non-negative kernel function K. The bandwidth h controls the bari-vari trade-off of the estimated Riemannian metric tensor and its derivatives.
Altering the link function y and the loss function Q in (3.7) enables flexible local regression estima- tion for models in Example 3.1.
Example 3.2. Consider the following loss functions for estimating the metric tensors and the Christoffel symbols when data drawn from model (3.2)-(3.4), respectively.
1. Continuous noisy response: use squared loss Q (y. (4) = (y - 11)² with g being the identity link function so JAu ?/--
2. Binary (dis)similarity response: use log- likelihood of Bernoulli random variable (v, (4) y log μ + (1 v) log (1 -m),
(3.11)
and g the logit link, so juu = logit (m)) ). The model becomes a local logistic regression.
3. Binary relative comparison response: apply the same loss function (3.11) and logit link as in the previous scenario, but here we formulate the linear predictor based on dist (Xuo; Xuz)² dist (Xao, X,1)² 1ha1 - 1/22 and
Nh = و (Ma 772)) Locally, the difference in squared distances approximated by
(3.12)
(60,0-10,0-1 5y -2 (
Tul Thu2
(3.13)
+ 0 (6 6112%1 - 8t 0 (8000 - 6
iske
for 8/ = X -pp and 8% 82,0-2 = 8iz2 - 5,00 i = I,..., d. Here Th1 and 7a2 are constructed in analogy to (3.6) using (X Xw1) and (X X22) 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 (Fan Gijbels, 1996), accommodating a wide variety of data. To efficiently estimate the metric on the entire manifold M, we apply procedure based discretization and post-smoothing, as detailed in Appendix B.
4 BIAS AND VARIANCE OF THE ESTIMATED METRIC TENSOR
This subsection provides asymptotic justification for model (3.2) with E (Yu|Xaao X.1) = dist (X.o, x under the squared loss Q(,,)) v) (μr y)² and the identity link g(H) μ.. The estimator we analyzed here fits local quadratic regression without intercept and approximates the squared distance by a simplified form of (3.6):
dist (X.o, X,1)² =s nh = 5 1000
(4.1)
for u = I,...,N. lation into matrix form. Denote the local design matrix and regression coefficients
N. Given a suitable order the indices i,j for vectorization, we rewrite the formu-
D. = (6%.0-18%00 1
8 8i 8y
gd 10.0-10 5% "", B B= -=
30
8
5

---

Under review as a conference paper at ICLR 2023
so that the linear predictor 7/u, B. Further, write D .D.....DN) DN)", y=
= =Y .Y Yv)² >
......
and W = diag (W1,.. WN), with weights Wn specified in (3.10). The objective function in (3.7) becomes (Y - (D3)) W (Y -y-De - and the minimizer is , (D"WD) D"WY, for which we will analyze the bias and variance.
To characterize the asymptotic bias and variance of the estimator, assume the following condi- tions are satisfied in neighborhood of the target p. These conditions are standard and analogous to those assumed in local regression setting (Fan & Gijbels, 1996).
(A1) The joint density of endpoints (Xxo, X.1) ) is positive and continuously differentiable.
(A2) The functions Gıj: 1 are C²-smooth for i, k = I,...,d.
(A3) The kernel K in weights (3.10) is symmetric, continuous, and has a bounded support.
(A4) sup, var (Xu, Xu1) cc.
Proposition 4.1. Under (A1)-(A4), bias (3|x) = Op (h²), var (aIx) = O, (N-1 as + 0 and 242+2d > oc, where X {(x,o, ) is the collection of observed endpoints.
The local approximation (4.1) is similar to local polynomial estimation of the second derivative of a varvare squared geodesic distance function, explaining the order of in the rates for bias and variance.
SIMULATION
We illustrate the proposed method using simulated data with different types of responses as de- scribed in Example 3.1. We study whether the proposed method well estimates Riemannian geo- metric quantities, including the metric tensor, geodesics, and Christoffel symbols, Additional details are included in Appendix C of the Supplementary Materials.
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 (a' (x',...,,"") rd)) as G.j =
(1 + E1+E xb -2 (1=)); for i,j 1,..., d. Under the additive model (3.2) in Example 3.1, we considered either noiseless or noisy responses by setting (p) = or (p) > 0 respectively.
Experiments were preformed with d 2 and the finding are summarized in Figure 5.1. For contin- uous responses, the left panel of Figure 5.1a visualizes the true and estimated metric tensors via cost ellipses (A.1) and the right panel shows the corresponding geodesics by solving the geodesic equa- tions (2.2) with true and estimated Christoffel symbols. The metrics and the geodesics were well estimated under the continuous response model without without additive errors, where the esti- mates overlap with the truth. Figure 5.1 evaluates the relative estimation errors G - G ||]|
and | | |
w.r.t. the Frobenius norm (A.2) for data from the continuous model (3.2).
For binary responses under model (3.3), Figure 5.1c visualizes the data where the background color illustrates A. Figure 5. and left panel of Figure 5. suggest that the intercept and the metric were reasonably estimated, while the geodesics are slightly away from the truth (Figure 5. la, right). This indicates that the binary model has higher complexity and less information is being provided by the binary response (see also Figure Cb in the Supplementary Materials).
5.2 RELATIVE COMPARISON ON THE DOUBLE SPIRALS
A set of x 104 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 (Xan- Xx1, Xx²) to generate relative comparison response Y, as follows. For u = I,...,N, Y, = 1 if Xwo, X,, belong to the same latent class and X no, X belong to different classes;

---

Under review as a conference paper at ICLR 2023
geadesics
0.6
2
0.3
0.0
-2
-0.3
-2
x1
-0.6
-2
-1
2
0.50
0.75
1.00
1.25
1.50
intercept
10 2030
bruth
estimated with additive
binary response
method
estimated with
estimated with binary response
(c) Simulated data, where line segments show pairs endpoints colored according to their bi- nary responses. The background shows the value of ft
(a) Estimated and true metric tensors using ellipses representation (left) and the geodesic curves (right) starting from (1,0) with unit initial velocities pointing to 1-12 o'clock directions.
Christoffel Symbol
estimated intercept
2
1
2
0
-2
2
relative
relative error
error
$200
080'0
$0000
00t'
SCU
o50
100
relative erro
$20
(b) Relative errors in term of Frobenius norm the estimated ten-
sors for the continuous response model (3.2) with additive error.
(d) Errors for estimating A with binary responses (3.3).
Figure 5.1: Simulation results for 2-dimensional sphere under stereographic projection coordinate.
otherwise Y O. Figure 5.2b shows portion of the N 6, 965, 312 comparison generated, where the hollow circles in the middle of each wedge correspond to X.v-
Here, contrast of the two latent classes induces the intrinsic distance, the distance is larger across the supports of the two classes and smaller within a 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 were 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 dis- tance and metric induced by the class membership discussed above. Meanwhile, the estimated geodesic curve unveils the hidden circular structure of the data support as shown in Figure 5.2d.
NEW YORK CITY TAXI TRIP DURATION
We study the geometry induced by taxi travel time in New York City (NYC) during weekday morn- ing rush hours. New York City Taxi and Limousine Commission (TLC) Trip Record Data was

---

Under review as a conference paper at ICLR 2023
(d) Geodesics starting from (1. with initial velocity pointing to 9 o clock under estimated metric.
response
latent class
class un A
class
(a) The points and their la- tent classes
(b) Relative comparison
(c) Estimated metric
Figure 5.2: Simulation results for relative comparison on double spirals. Gray curves (solid and dashed) in the background represent the approximate support of the two latent classes. In (b), the tiny circles are Xuo, each with two segments connecting to X,a and Xu2, colored according to Y,,
accessed on May 1st, 2022' to obtained business day morning taxi trip records including GPS coor- dinates for pickup/dropoff locations (X Xx1) and trip duration Y,. Estimation to the travel time metric was drawn under model (3.2) with Q 14) = (y na)² and g(p) μ..
Figure 6. la shows the estimated metric for taxi travel time. The background color shows the Frobe- nius 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 to 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 travel a unit distance along different directions. This suggests that in Manhattan, it takes longer to drive along the east-west di- rection (narrower streets) compared the north-south direction (wider avenues).
Geodesic curves in Figure 6. show where a 15-minutes taxi ride leads to starting from the Empire State Building. Each geodesic curve corresponds to one of 12 starting directions (1-12 o' o'clock). Note that we apply continuous Riemannian manifold approximation to the city, so the geodesic curves provide approximations to the shortest paths between locations and need not conform to 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 non-com con-constant curvature 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 big challenge apply nonparametric methods to data sources like images and audios. However, it often found that apparent high-dimensional data actually lie close some low-dimensional manifold, which is utilized by manifold learning literature produce reasonable low-dimensional coordinate representations. The proposed method can then be applied the resulting diweensinii coordinates as in the following MNIST example.
We embed images in MNIST to a 2-dimensional space via tSNE (Hinton & Roweis, 2002). Simi- larity 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 data induced by this similarity measures.
The geodesics estimated from our method tend to minimize the number of switches between labels. For example, the geodesic in Figure 7.1 remains "4" (Ist row of panel (b)) throughout, while the straight line on the tSNE chart translates to path of images switching between "4" and "9"
"https//www// wyy.gg//iiccee/aaoo//tttttrroooe
Data format changed after our download.
After rescaling, see Subsection C.4.
8

---

Under review as a conference paper at ICLR 2023
Morning rush hour traffic
Geodesics from the Empire State Building
40.85/M
MIULA
Frobenius nom
time (minutes))
40,80/N
15
0.5
0.3
40.79/N
0.2
40./O'N
74.00
73,95/W
(b) Geodesics correspond to 15 5--minute taxi rides from the Empire State Building heading to 1-12 o c'accc.
(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.
(2nd row of panel (b)); similar phenomenon occurs for geodesics B and C (corresponding to 4th and 7th rows in (b)). Also, our estimated geodesics produce reasonable transition and reside in the space of digits, while unrestricted optimal transport (3rd, 6th, and 9th rows of panel (b)) could produce unrecognizable intermediate images. Our estimated geodesic is faithful to the geometric interpretation that a geodesic is locally the shortest path.
4
9
9 4
4
4
3 3 3 3
3
3
3
3
3
3
KSNEZ
&
8
3
3
3
3
3
3
5 3 3 3
3
3
3
S 5 5 S 5
S
5 5 8 8 8
8
8
?
S 5 5 5
3 3
3 8
8
?
S S S
S
5
F 7
(b) Image transitions corresponding to A, B, and C (a). Ev- ery 3 rows correspond to a set of paths sharing the same pair of starting and ending images, where the first, second, and third rows correspond to the estimated geodesics, the straight lines the chart, and the optimal transport (path shown in (a)), respectively.
ISNE.1
(a) The estimated geodesic curves (solid red) and straight lines the chart (dashed black).
Figure 7.1: Geometry induced by a sum of Wasserstein distance and same-digit-ory -o--ooo-otott-t--- indicator.
DISCUSSION
We present a novel framework for inferring the data geometry based on pairwise similarity mea- sures. Our framework targets data lying on a now-dimensional manifold since observations need to be dense near the locations where we wish to estimate the metric. However, these assumptions are general enough for our method to be applied to manifold data with high ambient dimension in com- bination with manifold embedding tools. Context-specific interpretation of the geometrical notions, e.g., Riemannian metric and geodesics has been demonstrated in the taxi travel and MNIST digit examples. Our method has the potential to be applied many other application communities, such cognition and perception research, where psychometric similarity measures commonly made.

---

Under review as a conference paper at ICLR 2023
REFERENCES
Aurélien Bellet, Amaury Habrard, and Mare Sebban. Metric learning. Synthesis Lectures on Ar- tificial Intelligence and Machine Learning, 9():1-151, January 2015. ISSN 1939-4608. doi:
10.2200/S00626EDII VO1 2O1000001111111A1I Y201501AIM030.
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):3376, June 2004. ISSN 0148-9267. doi: 10. 10.116200148222422311222
Dong Chen and Hans-Georg Georg Müller. Nonlinear manifold representations for functional data. The Annals of Statistics, 40(1):1-29, February 2012. ISSN 0090-5364, 2168-8966. doi: 10.1214/
11-A0S936.
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.
Chopra, R. Hadsell, and LeCun. Learning a similarity metric discriminatively, with application to face verification. In 2005 IEEE Computer Society Conference Computer Vision and Pattern Recognition (CVPR' 005, volume 1, pp. 539-546 vol. 1, June 2005. doi: 10.1109/CVPR.200 202
Ronald R. Coifman and Stéphane Lafon. Diffusion maps. Applied and Computational Harmonic Analysis, 21(1):5-30, July 2006. ISSN 10635203. doi: a.001///.ccha220066 00.000004006....
Manfredo Perdigão do Carmo. Riemannian Geometry. Mathematics. Theory Applications. Birkhäuser, Boston, 1992. ISBN 978-0-8176-34900 778-3-7643-3490-11
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, 1st ed edition, 1996. ISBN 978-0-412-98321-4
Charles Fefferman, Sergei Ivanov, Matti Lassas, and Hariharan Narayanan. Reconstruction of Riemannian manifold from noisy intrinsic distances. SIAM Journal on Mathematics of Data Science, 273):770-808, September 2020. doi: 10.1137/1M112629X.
Andrea Frome, Yoram Singer. Fei Sha, and Jitendra Malik. Learning globally-consister local distance functions for shape-based image retrieval and classification. In 2007 IEEE 11th In- ternational Conference on Computer Vision, pp. 1-8, October 2007. 4408839.
doi: 10. 10910//CC2220077.
T. Hastie and R. Tibshirani. Discriminant adaptive nearest neighbor classification. IEEE Transac- fions on Pattern Analysis and Machine Intelligence, 18(6):607-616, June 1996. ISSN 01628828. doi: 10.1 1..11993..506411.
Søren Hauberg, Oren Freifeld, and Michael J. Black. A geometric take on metric learning. In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger (eds.), Advances in Neural Informa- fion Processing Systems 25, pp. 2024 2032. Curran Associates, Inc., 2012.
Geoffrey E Hinton and Sam Roweis. Stochastic Neighbor Embedding. Advances in Neural Information Processing Systems, volume 15, MIT Press, 2002.
Guido Kraemer, Markus Reichstein, and Miguel D. Mahecha. dimRed and coRanking Unifying Dimensionality Reduction in R. The R Journal, 10(1):342-358 2018. ISSN 2073-4859.
Serge Lang. Fundamentals of Differential Geometry. Graduate Texts in Mathematics. Springer- Verlag, New York, 1999. ISBN 978 1778-0337--88593-0 doi: 1-1197777---4---5111-
Tam Le and Marco Cuturi. Unsupervised Riemannian metric learning for histograms using Aitchi- son transformations. In Proceedings of the 32nd International Conference on Machine Learning, pp. 2002-2011. PMLR, June 2015.
G. Lebanon. Metric learning for text documents. IEEE Transactions on Pattern Analysis and Ma- chine Intelligence, 4844)897-500. April 2006. ISSN 0162-8828. doi:
10

---

Under review as a conference paper at ICLR 2023
Guy Lebanon. Learning Riemannian metries. In Proceedings the Nineteenth Conference on Un- certainty in Artificial Intelligence, UAI'03 03, pp. 362-369, San Francisco, CA, USA, 2002. Morgan Kaufmann Publishers Inc. ISBN 978-0-12-70566421
John M. Lee. Introduction Smooth Manifolds. Number 218 in Graduate Texts in Mathematics. Springer, New York; ; London, 2nd ed edition, 2013. ISBN 978-1-44199998111 978-1-4419-9982.
John M. Lee. Introduction to Riemannian Manifolds. Springer Berlin Heidelberg, New York, NY, 2018. ISBN 978-3-319-91754-2.
Clive Loader. Local Regression and Likelihood. Statistics and Computing. Springer, New York, 1999. ISBN 978-0-387-98775-0.
Francesca Mazzia, Jeff R. Cash, and Karline Soctaert. Solving boundary value problems in the open source software Package bvpSolve. Opuscula Mathematica, 34(2):387, 2014. ISSN 1232-9274. doi: 117449900pMaath220444422222
Wolfgang Meyer. Toponogov's theorem and applications. Lecture Notes, Trieste, 1989.
Dominique Perraul-Joncas and Marina Meila. Non-linear dimensionality reduction: Riemannian metric estimation and the problem of geometric discovery. arXiv:1305.7255 [stat], May 2013.
R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statis- tical Computing, Vienna, Austria, 2022.
Deva Ramanan and Simon Baker. Local distance functions:
taxonomy, new algorithms, and an evaluation. IEEE Transactions Pattern Analysis and Machine Intelligence, 33(4))794-806, April 2011. ISSN 1939-3539, doi: 10.1109/PPAM.221002222
Sam T. Roweis and Lawrence K. Saul. Nonlinear dimensionality reduction by locally linear embed- ding. Science, 29065500):2323-222222 December 2000. doi: 2..112//cciencee290555500
Dominic Schuhmacher, Björn Bähre, Carsten Gottschlich, Valentin Hartmann, Florian Heinemann, and Bernhard Schmitzer. transport: Computation of Optimal Transport Plans and Wasserstein Distances, 2022.
Matthew Schultz and Thorsten Joachims. Learning a distance metric from relative comparisons. In Advances in Neural Information Processing Systems, volume 16. MIT Press, 2003.
Karline Soetaert, Thomas Petzoldt, and R. Woodrow Setzer. Solving differential equations in R Package deSolve. Journal of Statistical Software, 33(9), 2010. ISSN 1548-7660. doi: 10.18637/ jss.v03 vss.v033.
Juan Luis Suárez, Salvador García, and Francisco Herrera. A tutorial on distance metric learning: Mathematical foundations, algorithms, experimental analysis, prospects and challenges. Neuro- computing, 425:300-33 322, February 2021. ISSN 0925-2312. doi: wwt.b.....coomm200.00.0.
Joshua B. Tenenbaum, Vin de Silva, and John C. Langford. A global geometrie framework for nonlinear dimensionality reduction. Science, 29065500):2319-2333, December 2000. doi: 10. 1126/science.220.550000
Kilian Q Weinherger and Lawrence K Saul. Distance metric learning for large margin nearest neigh- bor classification. Journal of machine learning research, 10(2), 2009.
Eric Xing, Michael Jordan, Stuart J Russell, and Andrew Ng. Distance metric learning with ap- plication to clustering with andeiinformation. In S. Becker, S. Thrun, and K. Obermayer (eds.), Advances in Neural Information Processing Systems, volume 15. MIT Press, 2002.
11

---

Under review as a conference paper at ICLR 2023
Supplementary Materials
CONTENTS
Introduction
2 Background in Riemannian Geometry
3 Local Regression for Similarity Measurements
3.1 Probabilistic Modeling for Similarity Measurements
3.2 Local Approximation of Squared Distances
Bias and Variance the Estimated Metric Tensor
5
Simulation
5.1
Unit Sphere
5.2 Relative Comparison the Double Spirals
6
New York City Taxi Trip Duration
7
High-dimeensional Data: An Example with MNIST
8 Discussion
A Additional Definition
12
Implementation Notes
13
C Additional Experiment Details
13
C.1 Unit Sphere and Round Metric
14
C.1.1 Bandwidth Selection
14
C.2 The Double Spirals
14
C.3 NYC Taxi Trips
16
C,4 The MNIST Example
17
Spread of Geodesies
17
D.I Proofs
19
Asymptotic of the Estimated Metric Tensor
22
A ADDITIONAL DEFINITION
A cost ellipse visualizes the metric ellipse
8, = (x xd)
Σ
(x*
(dd =p²
pd)
(A.1)
-==
12

---

Under review as a conference paper at ICLR 2023
for constant O, which shows, approximately, the intrinsic distance the manifold when traveling unit length on the coordinate chart along each direction. More precisely, shows the norm of tangent vectors v²0, E T..M subject to Σi=1 v'v² = r² at pointing to the corresponding direction. For example in d 2-dimensional manifold p O, under G diag (A1. X2) with Å A and Thus, the direction in which the ellipse is larger corresponds to the direction of the larger geodesic distance. One can see (A.1) as the "inverse" of the Tissot's indicatrix, where the latter shows local equidistance contour to the ellipse's center.
1. The long axis, (+vAi, 0), is the
of the tangent vector ==
Frobenius norm of tensors is denoted ||-1| defined as
1/2
1/2
).
|
Σ GiyGis
|r
Σ f,jk==
1
(A.2)
for metric tensor G and Christoffel symbol r.
B IMPLEMENTATION NOTES
(R Core Team, 2022) package is developed implement the proposed methods and nu- merical experiments.
We utilized an efficient procedure to obtain estimates Cig over the entire manifold M follows. We first obtain estimates (Pm( over dense grid of points p1 P2:........ . PN M by following (3.7)-3.9).. Next, the estimate Gzz(z) at an arbitrary x E M is obtained by the post-smoothing estimate
Gif(z) ENew K (1 - Pm h Gig(Pa)
INuig K (1|x Pxll /his)
for some kernel and bandwidth hps > O. We also use local regression (Loader, 1999) for post- smoothing. The grid for the examples (Section 5, Section 6, and Section 7) 128 x 128 for the unit sphere, and 80 x 80 for the double spirals, 250 x 250 meters cells for the New York taxi example, and 64 x 64 for the MNIST example.
The estimated geodesics are computed by numerically solving ordinary differential equations sys- tem, either given the start point and initial velocity, or given the start and the end points. It suffices to notice that the geodesic equations (2.2) are equivalently written as, after plugging-in the estimated Christoffel symbol I.,
(t) (t). (t) - (t)u² (t) rt()),
for i,j,h = I,...,d. Here, i&k (1) is the value of the estimated Christoffel symbol at point Further supplying initial condition 70(0) pb- v0) = v6. = 1..... d for point Pa E M and tangent vector uo E p. constitute an initial value problem, whose solution reflects the geodesic curve starting from Po with initial velocity On the other hand, supplying boundary condition 70) = ph- ((1
(7
i = l,...,d for Po, E M constitute boundary value problem, whose solution reflects the geodesic curve from Po Pl. we use deSolve (Soetaert et al., 2010) and bvpSolve (Mazzia et al., 2014) for initial value problems and boundary value problems respectively. See reference therein for further details of numeric solution to ODE.
C ADDITIONAL EXPERIMENT DETAILS
This section provides further detail to complete Section 5, Section 6, and Section the main text including how generated the simulated data, and more figures.
13

---

Under review as a conference paper at ICLR 2023
C.I UNIT SPHERE AND ROUND METRIC
The stereographic coordinate of the d-dimensional sphere St identifies points on the sphere by map- ping it to its stereographic projection in R from the north pole. The round metric on the sphere 8 is the metric induced by embedding of S + R For detail, see for example page 30 of Lee (2013) and chapter of Lee (2018). We generated endpoints X X,1 uniformly in the coordinate chart (-3. 3) x (-3, 3), then pair the endpoints so that the difference in coordinates of the endpoints 18 5
X X would not exceed 0.2 for 1, N.
More precisely, data were generated on d = 2-dimensional sphere under stereographic projection coordinate. total of = 5 x 100 pairs of endpoints with X X E (-3,3) were generated subject to XX 40 X < 0.2 for all u ,,...,:;; ratio, we set (p) = a ~ x 10 for all p. which is approximately one-tenth of the marginal expectation of squared distance, i.e., a Re E dist (Xwo, X.i)² /10.
1,...,d d. For a reasonable signal-to-noise
For simplicity of presentation, we scaled the distance for the binary similarity response model (3.3). More precisely, we use dist, (:-))
/c dist (-.-) induced by the scaled metric Gij. cGi for some constant and i,j L d. The experiment here used c = 300. Intuitively, the constant e regulates the signal-to-noise ratio without changing the form of geodesics. Given the endpoints, smaller leads to a smaller value of geodesic distance and hence smaller variation in the linear predictors !/u. so the response Y, will take less influence from the distance, representing a higher amount of noise. Then fo (p) was set to be the average local squared distances within local neighborhood p under the scaled distance.
In the end, the responses were generated following (3.2) and (3.3) respectively.
In addition, Figure C.I illustrates the relative Frobenius error for estimated tensors using noiseless or binary responses.
C.I.1 BANDWIDTH SELECTION
Like local regression, the proposed method relies neighborhood specification for opti- mal bias-variance trade-off. The simulation in Subsection 5.1 uses the rectangular kernel K(x) (x) for (3.10), where is the indicator function, the estimation only utilizes observations with endpoints Xan, Xu1 are both lying in the neighborhood Up (x''....."") :12 ppll h l... d] of the target point p.
train-test set scheme for driven bandwidth selection. To simplify computation, we only considered additive error under (3.2). A 16 x 16 grid pı...., P256 E (-3.3) x (-3.3) were used as target points where metric tensors observations that were within close proximity to the grid. Estimation of the tensors were computed w.r.t. bandwidth h utilizing a train set containing Numin 400158 (approximately 80% of the data)
We propose a
estimated, with test set of N 31246
randomly selected observations outside of the test set. For the test set, Yu,test fuu,es were then computed under identity link by plugging the estimated tensors into (3.6).
The bandwidth minimizing the squared loss EN (Yu.eatt - Yu, /Nesst is then chosen. The proposed bandwidth selection resulted in h = 0.18 according to the loss shown in the left panel of Figure C.2, which corresponds to an effectively local sample size of approximately 1000 in the training step (right panel of Figure C.2). These tuning parameter choices were applied in the results shown in Subsection 5.1. Other bandwidth selection methods developed for local regression (see e.g. Fan & Gijbels, 1996, section 4. 10) also be adopted here.
C.2 THE DOUBLE SPIRALS
Define class of spiral functions S R > R² with
t + (cos(5t + d), t sin + Φ)).
The underlying spirals for class A and B are SA S(, Φ = 0) and Sn S(, Φ TT) respectively. For endpoints, first generate points X,, independently and uniformly on SA or Sg. then the endpoints are generated following X,, 0.157m for i.i.d. standard Gaussian random variables
14

---

Under review as a conference paper at ICLR 2023
metric
Christoffel ymmbol
2
2
-2
-2
1
2
-2
1
2
X1
x1
value
value
(a) Noiseless
metric
Christoffel Symbol
2
1
X2
&
-2
-2
x1
x1
value
value
22'
OSO S2.
27
2.0
(b) Binary response
Figure C.I: Relative emrors w.r.t. Frobenius norm (A.2) of the estimated tensors with noiscless or
binary response for 2-dimensional sphere under stereographic projection coordinate chart.
7. m 1,..., 70000. Provided with those candidate endpoints, we pair them to form relative comparison subject to the restriction that X Xayl < 0.35 for i, = 1,2, 1,...,N. N. The responses Y, are then generated based on the class involving endpoints by their corresponding X on the spirals.
15

---

Under review as a conference paper at ICLR 2023
RMSE
average training size
0.0014
0.0013
4000
0.0012
3000
0.0011
2000
0.0010
1000
0.0009
0
0.1
0.2
0.4
0.2
0.3
bandwidth
Figure C.2: Root mean squared error (RMSE) for the test set (Left) and the average number of local training observations (Right) the bandwidth h varies.
For estimation, we used a larger local neighborhood Up = (<x,"²") |x p with / x 2/2 and weights w, 1 Xu for I,..., N to avoid degenerate esti- mates.
S h for i 1, 1.2)
Note that different starting points and initial velocities will generate different geodesics, not all resembling a spiral, as shown Figure C.3.
2
0
-2
label of geodesics
2 3
Figure C.3: Geodesics with different starting points and initial velocities under estimated metric, crosses indicate starting points.
C.3 NYC TAXI TRIPS
We focus on the 8,809,982 sensible records between a.m. to 10 a.m. on business days from May to September (summer months to hopefully avoid snow) of 2015 in New York City areas other than the Staten Island. Sensible in terms of GPS coordinates not falling in to the rivers, travel time not being several seconds, and that inferred traveling speed is not 120 mph, and e.t.c. We measure the cost to travel Y, by the squared trip duration (instead of the trip distance). For each target location p. estimation computed using trips among the M 10 closest pickup/dropoff endpoints
16

---

Under review as a conference paper at ICLR 2023
in the neighborhood Up { -p'1 kilometers for i 1, and weights given by 2.5 kilometers with the kernel K' being the density function of the standard normal distribution.
C.4 THE MNIST EXAMPLE
The dimension reduction is computed using R package dimRed (Kraemer et al., 2018). The Wasser- stein distance and optimal transport are computed using package r transport (Schuhmacher al., 2022). To show image transitions, weighted average is adopted to approximate the inverse of the tSNE embedding so as map the trajectories back image space, similar to, e.g., equation (3.9) of Chen & Müller 2),1),, but with Gaussian kernel and a sufficienlly small bandwidth.
To simplify computation, we only embed the first 10* images (half the entire data), and the resulting embedding coordinates were scaled (i.e., centered by the mean and divided by standard deviation). We generated N 10 comparison by selecting nearby points in the embedded space subject to Loxo Xw1 so I 0.75, whose response Y, computed based on the aime-digit-or- not indicator and 2 Wasserstein distance between the corresponding images:
dist (Xwo: x = Cdistwass (pic,o: pic,.) ) + (MM.))/ll.)>
where for u N,
Xuo, Xu1 E R are coordinates in the embedded space;
pic,o and pic,. are the 28 x 28 grey scale images;
1bl,0 and Ibl, are the image labels (0-9);
distwase (:--) is the 2-Wasserstein distance treating images 2-dimensional probability distributions;
1 (ovent) is the indicator for whether the event true (1) or false (0).
We multiply the Wasserstein distance by C = to balance the magnitude of the two summands, otherwise the later could be overly dominating.
Estimation was drawn under model (3.2) with squared loss (4) = (y- + and the identity link. We included the intercept term (8(0) in (3.6) to capture intrinsic variation. Since the dimensional reduction embedding map is necessarily similarity measures could share identical coordinates in the embedded space. Figure C.4a shows the estimated intercepts, which is larger among class boundaries, coherent greater variation in the similarity measure. Those few blank pixels indicate failure to obtain positive definite metric, which alleviated by averaging neighboring estimated values.
(y
injection, so that different images with non-zero
We also dropped the terms for Christoffel symbols from (3.6) for better numeric stability. Con- sequently, the estimated Christoffel symbols were computed by numeric differential following the definition. Results are similar if we include the Christoffel symbol terms in the linear predictor, but less stable
Figure C.4b shows the cost ellipses for addition visualization.
Notably, the proposal also work supplied with binary similarity measures using only the same-digit- or-not indicator (i.e., setting C' = 0 to remove the Wasserstein distance), and retains the "fewer label switching" tendency as illustrated Figure C.5. We see this as an real data example analogy to the double spirals (Subsection 5.2). We would also like to remark that not all geodesics are different from straight lines on the chart, and it is not guaranteed that geodesic must travel within the same class whenever possible, since its travel is jointly determined by the metric and where the endpoints are located.
D SPREAD OF GEODESICS
Here provide proof to the Proposition 3.1 in the main text, which characterizes the distance between geodesics departing from a same starting point. Proposition 3.1 is a result of combining Proposition D.I and Proposition D.2.
17

---

Under review as a conference paper at ICLR 2023
0
0
0
0
intercept
1.00
0
0
0.75
0:50
0.25
ISNE.1
SNEA
(b) The cost ellipses of estimated metric.
(a) Estimated intercept reflecting intrinsic local variation of the similarity measure.
Figure C.4: More figures for the induced geometry by adding Wasserstein distance and same-digit- or-not indicator.
C
7
7
9
7
7
94
O
O
O
O
O
O
O
O
2
2
2
2
2
2 2
1SNE.2
O
O
2
2
2
A
3 3 3 3 3
3
I
/
/
/
7
S
3 3
3
3
8
8
8
/
9 9 Z 2 2
2
2
2
3
0
0
0
b
9 9
6
(b) Image transitions (per row) corresponding to bows B,
C. and D in panel (a). Every 2 rows correspond to one pair of start and end images, where the first and second rows follow geodesies and the straight lines respectively.
-2
ISNE.1
(a) The geodesic curves (solid red) and straight lines (dashed black).
Figure C.5: Induced geodesics of only the same-digitoor-not indicator.
Proposition D.I (spread of geodesics). Let E M and u, w E TpMM be two tangent vector at p. Then the squared distance of separation satisfies Taylor expansion of
1 1a u) 3
dist (t)), (tw,) = t² ||0 - w|
as 0.
Here, R is the (1, curvvauuee tensor defined as
7 Vpx.?
18

---

Under review as a conference paper at ICLR 2023
where X,Y, Z are vector fields and [X,Y] XY - YX is the Lie bracket (c.f., e.g., Lee, 2018, page 385). Further, the Riemann curvature tensor is defined as
Rm(X, Y, W) W) (R(X, (
where W also a vector field. Note that R and Rm are both tensor fields, (R(r, w)w, v) (equivalently Rm((v, w, w, v)) are those evaluated at since w M.M. See Lee (2018), pp. 196- 199 for detail.
However, additional terms are introduced when computing via coordinate charts, as a result of ap- proximating the initial velocities v and w.
Proposition D.2 (approximation of velocity). p and 7(t) exp,(te) (tv) be the geodesic from p with initial velocity Given any local coordinate chart, write O. For 1,... d, denote 8² (t) (t) - 7 (0()) the difference in coordinate after traveling along 7, we have
any p E M, let v E T.MM be a tangent vector at
=t= +
where the remainder is
1
R (t)
8'"0 r
+
mnnTki 8T 8"88 (I I + a m + 0(+)
(D.1)
2t
6t
0, where I' and ar denote the the Christoffel symbols and their derivatives at p.
P
T,M
W
(tu,,(t)
spread of geodesics
(tw)
M
Figure D.I: visualization for the spread of geodesics as in Proposition 3.1. tangent space (blue planc) and tangent vectors are annotated in red.
D.I PROOFS
Proof of Proposition D.I. Similar results be found at Proposition 2.7 of do Carmo (1992), Proposition 5.4, of (Lang, 1999, IX, $5). We use the form presented by Meyer (1989). In the following, reproduce the proof the equation (9) of Meyer (1989) with some additional clarifi- cation.
Let (s() (sp) and (s))) exp, define family of curves
V(s,t) exp (texp- (o(s) (s)))
so that the curves V, : (is, t) are geodesics connecting (s) and (s) (c.f., e.g., proposiiin 5.19 and equation (10.2) of Lee (2018)), and that is a variation through geodesics V., Further, let T a,V, which is tangent field of velocities. Let E 0,V, which is a Jacobi field through geodesies V. that vanishes p. Denote H(s) dist (ro(s), (s)) || Is.t for any E [0. 1).
19

---

Under review as a conference paper at ICLR 2023
where the "Ia,t" means to take value point (ss, t). Then by the chain rules for covariant derivatives (see, e.g. Lee, 2018, chapter 4), we have
9,2 (DD,T, la,t d (T,T)
((())
ds
ds
d
2 H(s) (D + + latt
ds
d
8
H(s)
la,
ds
H(s) = (
(D,TT 1a.1.
ds
Note that Vo = for all t, so that -=, for all t, hence (0) = 0.
Note that V,, ( + (ss. (), + (ss, 0) and - V(s,1) 1) are geodesics; thus D,T 0 for all I, D,E
and D.Els.t ,,,=1
all 8. addition, by lemma 6.2 of Lee (2018), D,T = D,E.
By Jacobi equation, DE + R(E This means the vector field t + E' ==0, at is linear in t, together with E ==,,00 = v and E El = w, we can write
for all s, which implies DJ a=0 0 since 7 == 0.
E 0 = + t(w -
for E [0, 1]. Therefore D = D = w - v. which implies H" = 2 |/ w
Proceeding to the third order derivatives, observe that
H0) = 6 = 1==0,1.
and by proposition 7.5 of Lee (2018), D D.DEE D,D,E R(E,T)E, thus it suffices to show
D,E for all
(D.2)
in order to argue (0) 0. Since it is known that D, D,El==0,=0 0 == D, =00 it suffices to consider its derivative for (D.2). Use proposition 7.5 Lee (2018) repeatedly. we have
D,D,Els=a
D,D.D,F + Dr (R(T, E)E) -=
D.D +=00
-DD,,D,T- -E,,((EE)T 1=0,1
-D,,,DTT)==O,
D, (D,DIT R(E,T) I==0
D. (R(T, EET)) 0=0,,
where the last equation is due to D,T (c.f. e.g. proposition 4.15 of Lee (2018)).
0. Further, by chain rule of covariant derivative
D, (R(T, EYT) (VER) (T. EYT R(D,, R(T, 0, since a=, 0 for all t. Hence t + D.Els=o, is also linear
E)D,T,
which equals vector field, implying (D.2) and subsequently H0) 0.
For the fourth order derivative, note that (D.2) also implies that D,D,El.s =00, 0 and that D2T 0 for all f. Therefore,
((4) = ( 1--
Further,
D&T D.D,E
D. -D((D.D,E+RREET)))
20

---

Under review as a conference paper at ICLR 2023
= -D.,,,DDDDDDDEEEEEEEE.DE 1==0, Thus,
(D²T, l== (DD,D,D,E,,IY) + (.),, 1==0.--
Recall at 8 0, D2T 0.t D,Ells=ort
therefore
(R(E,D,) D, 1==0,t
Rom(E, D,T, E, D,T) 0-
Ron(v, w 0,0, w - v) + LRm(v, w - w u, w
Rom(e, w, w.
Further,
(D,D,D,E, D.T) _ (D,(D,DE,,D.T) (.T,E,,TTT) 1==0,-
(D,D,, D,T) 1s0,t
= D,D, (D,E, D.T) -= D3 (D,E, (², 1==0,0
where the second term in the last line vanishes since DJ Fl, = 0 and D2 El 0. Moreover, since the Levi-Civita connection torsion free, we have
(D.DD,EE, D.TT" la- = D,D, (D,E, D,T) I-= D.D, (D,E, -==0,-
which should be irrelevant to so that D& (D. E, D,T) l==0,t is linear in I. Yet
D3 ED, D,T) E,AE, D,T)) (D,E,
which vanishes at = 0 and = 0, 1. Hence D, (D, D,T) |== 0 for all f E [0. Combining those with the symmetries of Riemann curvature tensor leads the desired expansion.
Proof of Proposition D.2. Under the coordinate chart, we can write the geodesic curve as ( (t). (t)) for some smooth function 7) .'....??" Then for any i I,..., d, univariate Taylor expansion provides
7()) - 7000) -000 - -00
as 0, where 7',77 and 72 are the first, second, and third order derivative of 7t w.r.t. Note that the first derivative (0) v². and the geodesic equation and its derivative give
(0) = - " mn
T ) (0) = TT o (2T th MI
Plugging into the initial Taylor expansion gives the desired result.
Proof of Proposition 3.1 in the maintext. By Proposition D.1 and Proposition D.2, as t - 0, we have
2 E1u w|l
= (6 + tRi((t) - tt)i(t)) G.y (8 (o)(())))))))----))))( (s)(())
88-1% Gij + +208 (9)(() R((()) Gij + 0(1))
16-1% 1 + 8/ (8506 0/ (T + +(t4)
where
1
9 (1)
5mh T.
+
6t 0 828 (T
mn ki a ) O(t2)
2t
for a = 0, similar to (D.1). Note that 68 = (t) = (t), it suffices to keep only the first term in the ()), which is Ot)).
21

---

Under review as a conference paper at ICLR 2023
E ASYMPTOTIC OF THE ESTIMATED METRIC TENSOR
Now we discuss the variance and bias of the estimated metric tensors. For simplicity, use the squared loss Q(1, v) = (μ y) the identity link g(4) = 14, and exclude the intercept ((0) and the terms Pijk for derivative. Given a suitable order of the indices i, we rewrite (3.6) into matrix form. Denote
1,0-1 D, = (6,,0- 1,,0-10,0011----
8" 80 1,0-1
8n 10.0-18,,011)
B = (f 3 8())
dd
then the linear predictor 1/1 D B. Further, write
D= -D......, ,,,D,,""
Y= Y= -YI....,,, YN)", W diag (Wl,..,,N))
so the loss (3.7) becomes
(Y D W(Y - D33),
(E.1)
whose minimizer B WD) -1 WY. Therefore
bias (ÅID) (D" W',
var (öID) = (D"WWD) -1 D' DDD
where
r (E....)) - &f 1 Gig) IgnEN'
Σ diag ( var (YulX Xwo X ) ISnSN'
The assumptions (A1)-(A4) in the main text are reiterated here.
(A1) The joint density of endpoints X.o; X.1 is positive and continuously differentiable.
(A2) The functions Gij: Dt are C²-smooth for i,j. 1,...,d
(A3) The kernel K in weights (3.10) is symmetric, continuous, and has bounded support.
(A4) sup, var (Yu/X X.1) < 00.
Proposition F. Denote
SIN
Σ
who 0,0-1 500-15001100
5'2
S2N
Σ W var (Y,l Vulxuo, Xu1) 881
1210-18%000 00 4,0-19
84
SaN
=
W w
1/1 02,0-1100 Rus
where
Rg
Σ
(5h00 8 Gmr-
Under (A1), (A2), (A3), and (A4), and suppose that h - 0 and 24.dd - 0o, then
var SIN O (N 2-2 ESSNniisssss = =0 (Nh*),
F. = 0N 2d), var S2N = O (NA--6d), ESSN.isiniisi = (Nh'),
var SaN O (N 10-22
h and N 00. So
SiNdiigisi = O, (Nh*),
S2n =O (N/1-2d))
22

---

Under review as a conference paper at ICLR 2023
If further 2h oo, then
SaN = Op (Nh"),
h and N 00. Proof.
Write
Unihiiigh wuon 10,00-10000000 1 10,001%%0000
so
EU ugh fadgda
= =fn 2d II (K (82n/h) K (h))/h) 1
f(p' + ooo...., ,p + ot (41)d6 ... .o.
= / ($t0)) K (s (s sall (s (s)
(5 s(i))
(f(pt.... p²) + o(1) ds ds =0(1"),
where is the joint density of endpoints Xuo, Xu1, the second last equality is due to change of variables, and the last due to (Al). Similar argument implies
var Unshnépial 0 (h8-2d)
2d)
Ew EwgUUssitisk = O (h4-2d)
WyUlsstttti O (h8-od).
These rates apply uniformly over n, therefore by i.i.d. and that var Yul X21 is uniformly bounded,
ESIN, = 0 (NN"), var SIN = 0 (NN8
R.S.N.iiiiggi = 0
var S2N.ijiiiii = 0 (N18 B6))
Hence
SINdrisigi = ES Ni Op (V var SIN = O, (Nh"),
under h + 0 and 2h2dd 2d + cc. Similarly we have results for S2N
Next, write
Vusi = -Wodh.o Ra.
Note that
EV V
=f
=
2d
K (Cho/h) K (Bn1/h)) 15
Σ 5m (5% (5%
5th 6/.) Fxim
Iskdmmrd
f(p' + 81 610. ...
5 )d6 ... d.gd
155 K (s1u)
..... K (sa)) (s01
d1 sull (st0 <li)) x
SNO s (shosto s shi) FRim
IER.mrrd
0f (p) (s + + + o((1)
f(0) "...."))) + Σ
do de
Op
r=1
O(h))
23

---

Under review as a conference paper at ICLR 2023
where Fklm The Indeed, since the kernel K is symmetric, and the leading terms in the integrant is of fifth power of 8, thus with some abuse of notation, EV F ((s)s (f+ ho)))) ds 0(h)). Similarly
var Vasini 0(h
The rest of the results for S3N proceeds analogously to that of SINáiisiai and in
Under the conditions of Proposition E.I. Proposition E.2.
bias (g1x) = Op (h²). var (31x) =Op N7,
1
as h O and 212+2d - where X are the observed endpoints.
Proof. Note that i2ini are elements of DTDWD, where one pair of (2) index a row while one pair of (i3, i4) index column for i1, i2, is, I,..., d. Similarly S2Ndiizigg are elements of D'ED, and S3N elements of D Wr by Proposition 3.1, Applying Proposition E. leads to the result.
24


Reviews:
[{"review_id": "tFR4_2W29i", "reviewer": "Reviewer_gAuZ", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "The paper is concerned with learning a Riemannian metric such that geodesic distances matched (direct or indirect) measured distances. A general framework is presented that captures the common situations of (1) observations are distances, (2) observations are binary indicating if two points are similar or not, or (3) observations are triplets where one pair is more similar than others. It is then proposed to estimate the Riemannian metric through a local Taylor expansion. This is demonstrated to work well on fairly simple low-dimensional problems.", "summary_of_the_review": "The paper presents an interesting idea, but I find that the experiments are lacking.", "strength_and_weaknesses": "## Strengths\n* The derived Taylor expansion is $\\mathcal{O}(t^4)$ indicating that it is quite accurate as an estimate of a local metric.\n* The method seems both quite general and novel.\n\n## Weaknesses\n* I was unable to determine if the estimated Riemannian metric is always positive definite. It would be nice with a clarification on that.\n* As far as I can tell, Eq. 3.10 relies on Euclidean distances through the kernel functions. The main motivation of the work is that we should not trust the Euclidean distance, so if my reading is correct, then this seems a bit problematic.\n* I was unable to determine how geodesics were computed in practice.\n* The paper contains some simple toy experiments to verify that the method works and an application to travel times in Manhattan. While these experiments are nice, I miss a compelling application of the method. I believe that this line of work can be quite useful, but it would be good to demonstrate this to a greater degree.\n* I miss a more careful investigation in the synthetic experiments. For example, in the sphere example, it should be trivial to vary the sphere dimensionality and ask, \"does the method keep working?\" Given that only very low-dimensional problems are considered, I worry that the method does not scale. It would be nice with some evidence to counter such a belief.\n* The discussion section of the paper feels very rushed. Perhaps this is the usual time pressure of conference deadlines.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "M0h1noBpgWU", "reviewer": "Reviewer_LDTm", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "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 of the metric. 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. \n\nThe paper illustrates this solution on two simulated data (a unit sphere and a points on double spirals) and one real data set (time take by taxis between points in New York City. ", "summary_of_the_review": "\nThis is an interesting paper which solves a very specific problem. Estimating Riemannian metrics in general is much more ambitious and challenging but it assumes the knowledge of (noisy) pairwise geodesic distances which limits the scope and  makes the problem easier. ", "strength_and_weaknesses": "\nStrengths: \n\nSeems like a good idea to learn the local geometry of a manifold. In case that additional information (pairwise geodesic distances or labels or relative labels) is available and the points are relatively dense, then it seems reasonable to be able to learn the metric from the data. I don't know if any one has taken this approach in the past, but it seems quite novel to me. \n\nWeakness: \n•\tHow are the base points p (the basepoint for Eqn. 3.5) chosen in practice. The underlying approximation is valid only for a small t, so the base p should be close to the given data points. Does the method work if the given points are sparse?\n\n•\tWhat 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 I can’t think of too many other situations. \n\n•\tIt seems strange that the nonlinearity of the manifold is being learnt through a linear system of equations. \n\n•\tThe experiments are limited in scope. There are two simulation experiments – one on a unit sphere and one on spiral curves. The third example uses time of travel by taxis. A unit sphere is a special case where it is easy to learn the geometry through only the first order information (points and tangent vectors). It would be been nice to see an example of learning geometry of a more complicated manifold. \n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "HnJX-s3bq9s", "reviewer": "Reviewer_KGPf", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "main_review": "", "summary_of_the_paper": "To measure the difference between two samples $A$ and $B$,  this paper assumes the samples $A$ and $B$ are endpoints of two different geodesics starting from the same point $P$ in a small neighborhood and approximates their squared geodesic distances using Taylor extension which leads to a local linear regression model. Geometric characteristics like Riemannian metric and the Christoffel symbols serve as parameters of the proposed linear regression model. For the special case using squared loss, the authors established asymptotic bounds on the bias and variance of the estimator. Some intuitive results are shown through numerical experiments on unit spheres, spirals, and taxi trip duration. ", "summary_of_the_review": "I recommend \"reject, not good enough\" due to the scores regarding Clarity, Quality, Novelty And Reproducibility.\n[Quality and novelty of the proposed model, Score: 65/100]  [Quality regarding empirical use in modern CV&ML applications, Score: 20/100]\n[Clarity and Presentation, Score: 50/100] [Reproductivity, Score: 50/100] ", "strength_and_weaknesses": "[Strength]\n\n1). This paper proposes an interesting linear regression model which stems from a rigorous expression of Taylor approximation to the squared geodesic distances on a Riemannian manifold. \n\n2). Bounds on the bias and variance of the estimator are quantified in an asymptotic sense. \n\n[Weakness]\n\n1). Assumptions behind the local linear model may be too strict for real applications. \n\nThe proposed local regression model assumes the samples lie in small neighborhood of another point. However, in real applications using binary similarity measurements in computer vision, perception tasks, and recommendation systems, the compared samples may be far from each other and the assumption of localness may be violated. \n\n2). Insufficient empirical verifications for real applications.\n\nThrough experiments on unit spheres, spirals, and taxi duration data, the authors provide nice visualizations of the estimate geometric parameters. However, neither comparison results with other methods nor experimental results on modern CV&ML datasets are given. \n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject