Fast and Robust Shortest Paths on Manifolds Learned from Data

Georgios Arvanitidis†

Søren Hauberg†
Philipp Hennig‡
Michael Schober?
Technical University of Denmark, Lyngby, Denmark‡
University of Tübingen, Tübingen, Germany‡
Max Planck Institute for Intelligent Systems, Tübingen, Germany?
Bosch Center for Artificial Intelligence, Renningen, Germany†

Abstract

1.4
16

1.2
14

1

We propose a fast, simple and robust algorithm for computing shortest paths and distances on Riemannian manifolds learned from
data. This amounts to solving a system ofordinary differential equations (ODEs) subject to boundary conditions. Here standard
solvers perform poorly because they require
well-behaved Jacobians of the ODE, and usually, manifolds learned from data imply unstable and ill-conditioned Jacobians. Instead,
we propose a fixed-point iteration scheme forsolving the ODE that avoids Jacobians. Thisenhances the stability of the solver, while reduces the computational cost. In experiments
involving both Riemannian metric learningand deep generative models we demonstratesignificant improvements in speed and stability over both general-purpose state-of-the-art
solvers as well as over specialized solvers.

1

Introduction

A long-standing goal in machine learning is to build
models that are invariant to irrelevant transformationsof the data, as this can remove factors that are otherwise arbitrarily determined. For instance, in nonlinear
latent variable models, the latent variables are generally unidentifiable as the latent space is by design
not invariant to reparametrizations. Enforcing a Riemannian metric in the latent space that is invariant
to reparametrizations alleviate this identifiability issue,
which significantly boosts model performance and interpretability (Arvanitidis et al., 2018; Tosi et al., 2014).
Irrelevant transformations of the data can alternatively
Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics (AISTATS) 2019, Naha,
Okinawa, Japan. PMLR: Volume 89. Copyright 2019 by
the author(s).

12

0.8
10

0.6

8

0.4

6

0.2

4

0

2

-0.2
-1.5

-1

-0.5

0

0.5

1

1.5

Figure 1: A data manifold and the log-condition number of the Jacobian with fixed velocity (background).
High condition numbers can cause failures to converge
for off-the-shelf solvers.

be factored out by only modeling local behavior ofthe data; geometrically this can be viewed as having a
locally adaptive inner product structure, which can belearned from data (Hauberg et al., 2012). In both examples, the data is studied under a Riemannian metric,
so that can be seen as living on a Riemannian manifold.
While this geometric view comes with strong mathematical support, it is not commonly adopted. The
primary concern is that the computational overhead of
Riemannian models often outweigh the induced benefits. The main bottleneck herein are distance computations, which are often at the core of machine learningalgorithms. In Euclidean space, distance evaluations
require the computation of a norm. In kernel methods, distances are evaluated via the kernel trick, butmany interesting geometries cannot be captured with apositive definite kernel (Feragen et al., 2015). In bothcases, distances can be computed with negligible effort.
The distance between two points x, y on a Riemannian manifold is defined as the length of the shortest
path between x and y known as the geodesic. Computing the geodesic requires the solution of a boundary
value problem (BVP). These types of ordinary differential equations (ODEs) require specialized numericalmethods for their solution, as—unlike initial value prob-

Fast and Robust Shortest Paths on Manifolds Learned from Datalems (IVPs)—they cannot be solved via step-by-stepalgorithms like Runge–Kutta methods, and thus theirsolution is computationally taxing (Ascher et al., 1994).
Moreover, Riemannian manifolds learned from data
usually imply high curvature and an unstable Riemannian metric. The reason is that the metric is estimatedonly from finite data, so it changes irregularly fast.
As a consequence, the Jacobians of the ODE, which
are required in many off-the-shelf solvers, are often
ill-conditioned which causes additional problems andeffort (Ascher et al., 1994, §8.1.2). The example inFig. 1 shows that the Jacobian associated with a data
driven manifold is highly oscillatory. This implies thatstandard ODE solvers easily break down. Thus, specialized BVP solvers for shortest path computations
are required to make distance evaluations on manifolds
learned from data as fast and robust as in other models.
In this paper, we propose a novel method for computing shortest paths on Riemannian manifolds learned
from data. Combining the theory of smoothing
splines/Gaussian process regression (Wahba, 1990; Rasmussen and Williams, 2006) with fixed-point iterations(Mann, 1953; Johnson, 1972), we arrive at an algorithm
which is simple, fast and more reliable on challenging
geometries. This is achieved by not utilizing the Jacobian which is computationally costly and ill-behaved.
Our work is a significant improvement of an earlier
algorithm by Hennig and Hauberg (2014). Their algorithm is an early proof of concept for probabilistic
numerics (Hennig et al., 2015). We show below thatin the original form, it provably does not converge tothe true solution. However, their algorithmic structure is capable of converging to the true solution. We
demonstrate how their algorithm needs to be adopted
to improve solution quality and convergence speed, sothat the fixed-point algorithm can efficiently computethe shortest path on Riemannian manifolds.

2

A Brief Recap of Riemannian
Geometry

We start by defining Riemannian manifolds (do Carmo,
1992). These are well-studied metric spaces, where
the inner product is only locally defined and changes
smoothly throughout space.
Definition 1 A Riemannian manifold is a smoothmanifold M where each tangent space Tx M is equippedwith an inner product (Riemannian metric) ha, bix =|a M(x)b that changes smoothly across the manifold.
This inner product structure is sufficient for defining
the length of a smooth curve c : [0, 1] → M in the

Figure 2: A data manifold with a geodesic and its
tangent vector.
usual way as
Z 1p
Length(c) =ċ(t)| M(c(t))ċ(t)dt,(1)0where ċ = ∂t c denotes the curve velocity. The lengthof the shortest path connecting two points then constitutes the natural distance measure on M. The shortest curve is known as the geodesic, and can be foundthrough the Euler-Lagrange equations to satisfy a system of 2nd order ODEs (Arvanitidis et al., 2018),−1 | ∂vec[M(c(t))]−M(c(t))c̈(t) =ċ(t)2(ID ⊗ ċ(t) )
2
∂c(t)
|

∂vec[M(c(t))](ċ(t) ⊗ ċ(t)) , (2)−∂c(t)
where vec[·] stacks the columns of a matrix into a vectorand ⊗ is the Kronecker product.
Most numerical calculations on Riemannian manifolds
are performed in local tangent spaces as these are
Euclidean. Key operations are therefore mappings backand forth between the manifold and its tangent spaces.
A point y ∈ M can be mapped to the tangent spaceat x ∈ M by computing the shortest connecting curveand evaluating its velocity v at x. This is a tangent
vector at x with the property that its length equals
the length of the shortest path. By the Picard-Lindelöf
theorem (Picard, 1890), this process can be reversedby solving Eq. 2 with initial conditions c(0) = x andċ(0) = v. Mapping from the manifold to a tangent
space, thus, requires solving a boundary value problem,
while the inverse mapping is an initial value problem.
Practically, the BVPs dominate the computational
budget of numerical calculations on manifolds.
In the context of this paper we will focus upon
Riemannian manifolds which are learned from the data,and capture its underlying geometric structure. Thus,
we will consider for the smooth manifold the Euclidean
space as M = RD , and learn a Riemannian metric
M : RD → RD×D . To be clear, this simply changes
the way we measure distances, while respecting thestructure of the data. An illustrative example can beseen in Fig. 2.

Georgios Arvanitidis† , Søren Hauberg† , Philipp Hennig‡ , Michael Schober?

3

A Fast Fixed-Point Method for
Shortest Paths

In order to apply Riemannian models to interestingdata sets, we require a fast and robust method to solve
the boundary value problemc̈(t) = f (c(t), ċ(t)), c(0) = x, c(1) = y(3)

where f (c(t), ċ(t)) is the right-hand side of Eq. (2) andx, y ∈ M. Numerical BVP solvers typically replace
the analytic solution c(t) by an approximant c(t) whichis required to fulfill the ODE c̈(tn ) = f (c(tn ), ċ(tn )) on
a discrete mesh ∆ = {t0 = 0, t1 , . . . , tN −1 = 1} ⊂ [0, 1]of evaluation knots tn . Together with the boundary
conditions (BC) c(0) = x, c(1) = y, this results in a
D(N + 2)-dimensional nonlinear equation which can besolved for rich enough D(N +2)-dimensional parametricmodels. If the approximant is represented by the posterior mean µ(t) of a GP regressor GP(c̃(t); µ(t), k(t, s)),then it can be proven that a solution with small approximation error exists (Wendland, 2004, §11) undersuitable conditions on the kernel (Micchelli et al., 2006;
Rasmussen and Williams, 2006). This means we have tofind a (artificial) data set D with (N +2) D-dimensional
data points such that the posterior mean fulfills the
BVP on the discretization mesh ∆ and the BC on the
curve boundary B = {0, 1}. The data set D should bethought of as a parametrization for the solution curve.
One way to generate a solution is to define the auxiliary
function F (c) = c̈ − f (c, ċ) and apply a variant of the
Newton-Raphson method (Deuflhard, 2011) to find aroot of F . For example, this type of algorithm is also
at the core of Matlab’s bvp5c, a state-of-the-art BVP
solver. The convergence of these algorithms depends
on the Jacobians ∇c f, ∇ċ f of f at the evaluation knots
tn . In particular, Jacobians with big condition number
may cause Newton’s method to fail to converge if noprecautions are taken (Ascher et al., 1994, §8.1.2). On
manifolds learned from data, this is a common problem.
Furthermore, in practice these Jacobians are computedwith a computationally taxing finite difference scheme.
Thus, a method not based on Newton’s method should
be more suitable for the computation of shortest paths.
3.1

Method Description

As mentioned above, we model the approximate solution c(t) with the posterior mean µ(t) of a (multioutput) Gaussian processGP(c̃(t); µ(t), V ⊗ k(t, s))(4)

with (spatial) kernel k and inter-dimensional covariancematrix V ∈ RD×D . If the kernel k is sufficiently partially differentiable, this implies a covariance between

derivatives of c̃ as well (Rasmussen and Williams, 2006,§ 4.1.1), in particular

 m
dn
∂ m+n
d
c̃
(t),
c̃
(s)
=
V
k(t, s) (5)
cov
i
j
ij
dtm
dsn
∂tm ∂snfor the covariance between output dimensions i andj, derivatives m and n and spatial inputs t and s.
This (prior) model class is the same as in Hennig andHauberg (2014).
The boundary equations fix two parameters (0, x), (1, y)of the parametrization. The remaining N parameters(tn , z̈n ) approximate the accelerations c̈(tn ) of the true
solution c(tn ) at knot tn , i.e., z̈n ≈ c̈(tn ). The z̈n
are updated iteratively and we denote values at thei-th iteration with the superscript (i), e.g., c(i) (t) for(i)
the i-th approximation, z̈n for the i-th value of theparameter z̈n and so forth.
At iteration i, the approximation c(i) (t) is the predictive posterior GP
P (c̃(i) (t)) = GP(c̃(i) (t); µ(i) (t), k(i) (t, s))#"∂2k(B, B)
2 k(B, ∆)
∂s
G=V⊗
∂2
∂4
∂t2 k(∆, B)∂t2 ∂s2 k(∆, ∆)+ diag(0, 0, Σ, . . . , Σ)
(6)

h
i
∂2
ω | = V ⊗ k(t, B) ∂s
G−1
2 k(t, ∆)| 
x − m(0)µ(i) (t) = m(t) + ω | vec  y − m(1)  (i)z̈∆ − m̈(∆)
k(B, s)(i)|k (t, s) = V ⊗ k(t, s) − ω V ⊗ ∂ 2,∂t2 (∆, s)
and similarly for the velocity ċ(t) by forming the Gand ω accordingly. In Eq. (6), z̈∆ ∈ RN ×D represents
the accelerations, and k(B, ∆) ∈ R2×N is the matrix ofkernel evaluations at boundary points and evaluation
knots and similar for k(B, B), k(∆, ∆) and k(∆, B). Finally, Σ = εID is the identity matrix times a small(i)regularization parameter ε ≈ 10−10 , so µ̈(i) (tn ) → z̈n .
The rationale for its inclusion will be postponed to theend of Sect. 3.2 as it will become more apparent incontrast to the model of Hennig and Hauberg (2014).
Details for the components k(·, ·), m, ∆ in Appendix A.
Now, we proceed with the description of the algorithm.
Just like a root of the function F (c) = c̈ − f (c, ċ)solves the ODE, so does a fixed point of the mappingµ̈(i+1) (t) = f (µ(i) (t), µ̇(i) (t)). In particular, we canevaluate this mapping on the discretization mesh ∆ to
(i)
(i+1)
map z̈∆ to z̈∆ . The big advantage of this combination of parametrization and update scheme is the simplicity of obtaining closed-form iteration updates (Rasmussen and Williams, 2006, § 9.4). The vector field f is

Fast and Robust Shortest Paths on Manifolds Learned from Data

Step 1
3

15
1.4

10

2

1

1

velocity

acceleration

1.2

5

0.8

0

c2

0.6

-5

0

-1

0.4

-10

-2
0.2

0

-15

-3

-0.2

-20

-4
0

0.1

0.2

0.3

0.4

0.5

0.6

time

0.7

0.8

0.9

-1

1

-0.8

-0.6

-0.4

-0.2

Step 2

0

c1

0.2

0.4

0.6

0.8

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

time

(i)

Figure 3: Step 1 : With the current estimates of z̈∆ ( ) we generate using the GP model Eq. 6, the current
posterior curve µ(i) ( ) and velocities µ̇(i) ( ). Step 2 : Then, using the proposed fixed-point update scheme(i+1)
Eq. 7, we get the updated parameters z̈∆ ( ) . The algorithm iterates until k− f ( , )k small enough.
evaluated using the current iteration (µ(i) , µ̇(i) ) to yield
(i+1)
(i+1)
= f (µ(i) (tn ), µ̇(i) (tn )), and z̈nz̈n≈ µ̈(i+1) (tn )(i+1)
because ε → 0. Forming µ(i+1) , µ̇(i+1) from z̈n
only
requires two matrix-vector products (see Eq. (6)). The
process is depicted in Fig. 3.

is not stated. Interpreting the iteration as a fixed pointsearch is the key insight of this paper.
We suggest to apply a Mann iteration (Mann, 1953;
Johnson, 1972) process for the solution of (2) given byz̈∗n = f (µ(i) (tn ), µ̇(i) (tn ))z̈(i+1)= αi z̈∗n + (1 − αi )z̈(i)nnAlgorithm 1 The proposed fixed-point method.
Require: BVP f (c(t), ċ(t)), hyper-parameters ∆, ε,tolerance τ(i)
1: # Compute µ(i) (t), µ̇(i) (t) using Eq. 6 and z̈∆ .
(i)

(i)

2: Define: en , z̈n − f (µ(i) (tn ), µ̇(i) (tn ))2(0)

3: z̈n ← 0, n = 0, . . . , N − 1
4: i ← 0(i)
5: while ∃n : en > τ do
6:z̈∗n ← f (µ(i) (tn ), µ̇(i) (tn )), n = 0, . . . , N − 1
7:
for j = 0, . . . , 3 do
8:αj = 3−j
(∗,j)

(i)

z̈∆ ← αj z̈∗∆ + (1 − αj )z̈∆P (∗,j) P (i)10:
if n en ≤ n en then
11:
break
12:
end if13:
end for
(i+1)
(∗,j)
14:
z̈∆
← z̈∆
15:
i←i+1
16: end while
17: return GP(c̃(t); µ(i) , k(i) )9:Variants of this scheme have been repeatedly applied forthe creation of probabilistic differential equation solvers(Hennig and Hauberg, 2014; Chkrebtii et al., 2016;
Schober et al., 2014; Cockayne et al., 2016; Kersting andHennig, 2016; Teymur et al., 2016; Schober et al., 2018;
Kersting et al., 2018). Of these papers, only Hennig andHauberg (2014) points out that this can be updated
multiple times, but even there the connection between afixed point of the mapping and an approximate solution(7)

with “step sizes” αi ∈ [0, 1]. The results of Mann (1953);
Johnson (1972) only apply if αi = (i+1)−1 , however wefound a backtracking scheme to be effective in practice.
Algorithm 1 presents our method in pseudo-code where(∗,j)
z̈∆ denotes the tentative parametrization. Note how
the backtracking line search for αi (Lines 7-13) requires
half of the description.
Our method is similar to a recently proposed method
by Bello et al. (2017) that is based on the variational
iterative method 1 by He (2000). Bello et al. (2017)
proposed to use this scheme symbolically requiring acomputer-algebra system for its execution, which makes
it inapplicable to practical tasks. More details to these
related works can be found in Jafari (2014); Khuri andSayfy (2014).
3.2

Comparison with Hennig and Hauberg(2014)The proposed method is inspired by the previous workof Hennig and Hauberg (2014) and we make a directcomparison here.
The algorithm of Hennig and Hauberg (2014) is a
proof-of-concept probabilistic numerical method (Hennig et al., 2015) for solving boundary value problems. It
is structurally similar to other early probabilistic IVP
solvers of Chkrebtii et al. (2016) and Skilling (1991).
Since the publication of these early works, the field
1

calculus of variations not variational inference.

Georgios Arvanitidis† , Søren Hauberg† , Philipp Hennig‡ , Michael Schober?

rameters θ might require an exponential bigger
mesh (Vaart and Zanten, 2011, Thm. 10), theproperty of universality holds regardless of thehyper-parameter θ used to find the approximation. As a consequence, tuning hyper-parameters
is purely optional and certainly does not require
restarts. This also improves runtime significantly
as the Gram matrix needs only be to inverted once.
In practice, we have observed negligible solution
improvements after type-II maximum likelihood
optimization after the end of the algorithm.

has matured significantly, providing algorithms withnovel functionality (Hauberg et al., 2015; Mahsereciand Hennig, 2015; Oates et al., 2017; Xi et al., 2018)
and rigorous analysis (Briol et al., 2015; Chkrebtii et al.,2016; Schober et al., 2018; Kersting et al., 2018).
Their main idea is to treat the vector-field evaluations z̈n as noisy observations of the true, but unknown, second derivative c̈(tn ). For a concrete suggestion, they propose a Gaussian likelihood P (z̈n ) =N (z̈n ; c̈(tn ), Λn ). Together with a GP prior on c(tn ),they arrive at an inference algorithm. Heuristically,
they propose to add mesh observations sequentially,refine them iteratively for a fixed number of steps, andthey repeat the overall process until they find a set of
hyper-parameters of the GP which maximizes the data
likelihood of the final approximation.
However, the algorithm of Hennig and Hauberg(2014) cannot converge to the true solution in general. A convergent method is required to satisfyf (µ(i) (tn ), µ̇(i) (tn )) → µ̈(i) (tn ) as i → ∞. However,as Σn 6= 0 in Hennig and Hauberg (2014), µ̈(i) (tn ) is(i)
not an interpolant of z̈n (Kimeldorf and Wahba, 1970,
Thm. 3.2; Kanagawa et al., 2018, Prop. 3.6) implyingthat the true accelerations c̈(tn ) cannot be a parameterization in the model of Hennig and Hauberg (2014)
contradicting the fixed point requirement.
The same criticism could be applied to our model,
as we propose Σn = Σ = εID 6= 0. We have experimented with annealing schemes for this hyperparameter Σ(i) = i−1 Σ, but the benefit of ε > 0 for thestability of the Gram matrix G is bigger than induced
numerical inprecision, in particular when the tolerance
τ is considerably larger than ε.
Although the algorithm of Hennig and Hauberg (2014)
cannot converge, three insights and resulting modifications lead to our proposed method:
1. Hennig and Hauberg (2014) did not propose aprincipled scheme to determine the number of refinements S, but treated it as a hyper-parameterthat must be provided by the user. However,it can be easily checked whether the posteriormean µ(i) (t) fulfills the differential equation at
any point t. This not only removes a hyperparameter, but can also gives more confidence
in the returned solution. In principle, the er(i)ror en = µ̈(i) (tn ) − f (µ(i) (tn ), µ̇(i) (tn )) could
even be used to construct adaptive meshes ∆(i)
(Mazzia and Trigiante, 2004).
2. Using an universal kernel (Micchelli et al., 2006)
and a fine enough mesh ∆∗ , it is known that any
curve can be fitted. While sub-optimal kernel pa-

3. Currently, there is no analysis when or if
coordinate-wise updates offer improvements over(i)simultaneous updates for all parameters z̈∆ .
There is, however, a strong argument for updating
simultaneously: runtime. All predictive posterior
parameters can be pre-computed and kept fixed
throughout the runtime of the algorithm, if themesh is not adapted throughout. In particular, theregression weights ω | can be kept fixed, so each
update requires only two vector-vector products.
Finally, while our derivation does not make use of its
probabilistic interpretation, further steps in this direction could potentially unlock novel functionality which
has repeatedly been the case with other probabilistic
numerical methods (Xi et al., 2018; Oates et al., 2017;
Hauberg et al., 2015).

4

Experiments

In this section we demonstrate the advantages of our
method compared to Matlab’s bvp5c, and the algorithm in Hennig and Hauberg (2014) denoted in theexperiments as H&H. Since all the methods depend
on a set of parameters, we will come up with a default setting. For the proposed method, we will use
for the ∆ a uniform grid of N = 10 points including the boundaries. The corresponding noise term of
the points z̈ will be kept fixed to Σ = 10−7 ID . For
the GP we will use the Squared Exponential kernel2k(t, t0 ) = exp(−(2λ2 )−1 |t − t0 | ). We fix the amplitude V in a Bayesian fashion as Hennig and Hauberg(2014), and the length-scale λ2 ≈ 2−1 |tn+1 − tn | which
provides enough degrees-of-freedom while covering theentire interval at the same time. The prior mean is set
to the straight line m(t) = c(0) + t · (c(1) − c(0)), andthe derivative accordingly. For the method of Hennig
and Hauberg (2014), we use the same parameters as
for our method. We set the bounds on the Jacobian to
U = U̇ = 10 and run the method with one refinement
iteration. We set the maximum mesh size for the bvp5c
to 1000, and use for the starting mesh uniformly 10
points on the straight line connecting the boundary

Fast and Robust Shortest Paths on Manifolds Learned from Data
10 1

1.5
Ours Vs H&H
Ours Vs bvp5c

Curve Length

0.5

Elapsed time (sec)

2.5
1

2
1.5
1

10 0

10 -1

10

-2

0.5

0

0
-1

-0.5

0

0.5

1

0

0.5

(a) Data

1

1.5

10 -3

0

50

100

Curve lengths (ours)

Point index(b) Curve Lengths(c) Runtimes

150

200

Figure 4: Left: Generated data on a semi-circle, together with some challenging geodesics computed by our
method, and a point y. Middle: The curve lengths of the geodesics between the given data and the point y, on
the horizontal axis we have the point index, and the vertical lines represent the failures of the bvp5c to converge.
Right: The runtime for the corresponding geodesic problems.

4.1

Experiments with a Non-parametricRiemannian Metric

We first consider the case of Riemannian metric learning
as proposed by Arvanitidis et al. (2016). This canbe seen as a way to capture the local density of thedata, and thus to uncover the underlying geometricstructure. The metric at a given point is computedin three steps: 1) use a kernel to assign weight togiven data, 2) compute the local diagonal covariancematrix, and 3) use its inverse as the metric tensor.
More formally, the metric is
!−1
N
X
2
Mdd (x) =wn (x)(xnd − xd ) + ρ,(8)curves. Also, we see that the runtime of our method(Fig. 4c) is increased only slightly for the difficult
problems, while bvp5c is always slower and especiallyfor the difficult problems increases the mesh size to the
maximum and fails to converge. The performance ofthe H&H remains almost constant, since it is essentiallybased on the converge of the model’s parameters andnot the difficulty of the problem.
1
Fraction of failed geodesics

points. For all the methods we consider the resulting
2
curve as correct if kc̈(tn ) − f (c(tn ), ċ(tn ))k2 ≤ 0.1, ∀n.

0.8

Ours
H&H
bvp5c

Curly
0.6
0.4

Two-moons
0.2

n=1
2

kxn − xk2
where wn (x) = exp −
2
2σM

!

0

.

The parameter σM ∈ R controls the curvature of the
manifold, since it regulates the rate of the metric change
i.e. when the σM is small then the metric changes fast
so the curvature increases. The parameter ρ ∈ R>0 isfixed such that to prevent zero diagonal elements.
We generated 200 data points along the upper semicircle and added Gaussian noise N (0, 0.01) as it isshown in Fig. 4a. The ρ = 0.01 is kept fixed in all
the experiments. We set the parameter σM = 0.15,
and after fixing the point c(0) = y we compute thegeodesics to the given data. From the results (Fig. 4b)
we see that all three methods perform well when the
distance from the starting point is small. However,when distances increase, we see that only our methodmanages to find the correct curve, while bvp5c is notable to solve the problem, and H&H finds too long

Curly

Two-moons

Sphere-3D

SphereFigure 5: Failed geodesics.
Next, we generated three challenging datasets consisting of 400 points each. For the first one we generatea circle and we flip the lower half along the y axis, werefer to it as Curly in the results. For the second wemove the lower half of the circle such that to get thetwo moons. The third one is a 2-dimensional spherein R3 . Finally, we add Gaussian noise N (0, 0.01), andwe standardize to zero mean and unit variance each dimension. We keep the same parameters for the metric.
Note that the resulting manifold implies high curvature, so we increased the flexibility of the methods. For
the proposed model and H&H we used a grid of 50
points, and the maximum mesh size of bvp5c was set
to 5000. We pick randomly 40 points for each dataset,and compute the pairwise distances. In Fig. 5 we seeGeorgios Arvanitidis† , Søren Hauberg† , Philipp Hennig‡ , Michael Schober?

that our method manages to solve almost all of theshortest path problems, while bvp5c fails in almost halfof them. H&H is expected to fail in many cases, since
it is not designed to converge to the correct solution.

Mean elapsed time (sec)

40
Ours
bvp5c
H&H

30

20

4.2

10

0

2

5

10

20

Dimensions

D
Ours
bvp5c
H&H

2
0%
9%
69%

5
0%
24%
95%

10
0%
42%
100%

20
9%
50%
100%

Figure 6: Scalability in higher dimensions and failures.
Furthermore, we tested the scalability of the methods
with respect to dimensionality. We generate 1000 points
on a semi-circle in 2 dimensions, and standardized tozero mean and unit variance. We fix a point y and asubset S of 100 points. Then, for every D = [2, 5, 10, 20]
we construct an orthogonal basis to project the data
in RD , where we standardize the data again and add
Gaussian noise N (0, 0.01). Then, for each dimension D
we compute the geodesics between the y and the subset
S of the points. Keep in mind that the parameter σM =0.25 is kept fixed, so as the dimensions increase thesparsity of the data increase, and so does the curvatureof the manifold. In Fig. 6 we show the average runtime
for every dimensionality and the percentage of failuresfor each method. Our method remains fast in higher
dimensions, even if the curvature of the manifold isincreased. On the other hand, we observe that bvp5c
fails more often, which causes the overhead in theruntime. H&H remains fast, but, the resulting curves
cannot be trusted as the criterion of correct solution isalmost never met.

Mean data log-likelihood

-3
Ours
H&H
bvp5c

-3.5

-4

-4.5

-5

Also, we fitted a mixture of LANDs (Arvanitidis et al.,2016) using the three models, on the two moons dataset
generated by flipping and translating the data in Fig. 4a.
Note that we fix σM = 0.1 since we want our metric tocapture precisely the underlying structure of the data,
which implies that the curvature is increased. From theresults in Fig. 7 we see that the proposed solver fasterachieves higher log-likelihood. Additionally, in thesame time interval, it manages to run more iterations(dots in the figure).

0

1

2

3

4

5

Elapsed time (hours)

Figure 7: LAND experiment.

Experiments with a Riemannian Metricof Deep Generative Models

The Variational Auto-Encoder (VAE) (Kingma andWelling, 2014; Rezende et al., 2014), provides a systematic way to learn a low dimensional latent representation of the data, together with a generator thatlearns to interpolate the data manifold in the input
space. Usually, deep neural networks are used to model
the generator. These flexible models are able to compensate for any reparametrization of the latent space,which renders the latent space unidendifiable. Recently,
Arvanitidis et al. (2018) defined a Riemannian metricin the latent space, which is induced by the generator
and is invariant to the parametrization of the latent
space. This resolves the identifiability issue and makes
computations in the latent space parametrization invariant. More specifically, the VAEs utilizes a stochastic
generator that maps a point x from the latent space
X to a point y in the input space Y, and it consistsof two parts: the mean and the variance function asy(x) = µ(x) + σ(x) , where  ∼ N (0, Idim(Y) ) andis the pointwise multiplication. This stochastic mapping introduces a random Riemannian metric in thelatent space. However, as it is shown (Tosi et al., 2014)
we are able to use the expectation of the metric which
has the appealing form||M(x) = Jµ (x) Jµ (x) + Jσ (x) Jσ (x),(9)

where J stands for the Jacobian of the corresponding
functions. The interpretation of the metric is relatively
simple. It represents the distortions due to the meanfunction and the uncertainty of the generator.
In this context, for the data in the input space we generated the upper half of a 2-dim sphere in R3 , added
Gaussian noise N (0, 0.01) and scaled the data in theinterval [−1, 1]3 . Then, we trained a VAE, using forthe generator a simple deep network consisted of two
hidden layers with 16 units per layer, the softplus
as activation functions, and tanh for the output layer.
We used a 2-dimensional latent space, and the encoded
data can be seen in Fig. 8a. There, we show the computed geodesic for the 3 methods. Interestingly, we seethat the 3 resulting curves differ, and the estimatedFast and Robust Shortest Paths on Manifolds Learned from Data
10 2

3.5

Ours Vs H&H
Ours Vs bvp5c

Elapsed time (sec)Curve Length

3
2.5
2
1.5
1
Latent codes
Ours
bvp5c
H&H

Ours
bvp5c
H&H

10 1

10 0

10

-1

0.5
0

0

(a) Latent space

1

2

3

10 -2

0

200

400

600

800

Curve lengths (ours)

Point index(b) Curve Lengths(c) Runtimes

1000

1200

Figure 8: Left: The latent space together with the computed geodesics between two points. Middle: The curve
lengths, on the x axis we show the curve length of the proposed model. The results are sorted with respect to ourmodel. Right: The corresponding runtimes.
curve lengths are: proposed (2.52), bvp5c (3.65) andH&H (3.30). Our model, manages to find the shortest
path, which is a particularly curved path but the second derivative remains relatively smooth, while bvp5c
finds a simpler curve with larger length. This is not
surprising since bvp5c prefers solutions where the curveis smoother, while our method prefers curves where
the second derivative is smoother. In order to further
analyse this behavior, we randomly pick 50 points andcompute all pairwise distances. The results in Fig. 8bshows that the proposed method manages to find always the shortest path, while the other methods when
the distances increase, provide a suboptimal solution.
Comparing the runtimes (Fig. 8c) we see that our
method is faster in the simple problems, and has onlya small overhead in the difficult problems, however, it
manages always to find the shortest path.

Elapsed time (sec)

10 2
Ours
bvp5c
H&H

10 1

10 0

10

-1

10

-2

0

20

40

60

80

100

Point index

Figure 9: Runtime comparison
As a last experiment, we generated a 2-dimensional
sphere in R3 and moved the upper half by 1, and again
added Gaussian noise and scaled the data to [−1, 1]3 .
Instead of softplus in the hidden layers, we used thetanh activation functions which increases the curvature.
Here the latent space is 2-dimensional. We fix randomly
a point in the latent space and compute the geodesic to
100 randomly chosen points. As we observe from theresults in Fig. 9, as long as we compute the distance
between points of the same semi-sphere the runtimesof our method and bvp5c are comparable. However,when the points belong in different semi-spheres theruntimes increase significantly. The reason is thatcurvature increases dramatically when we cross partsof the latent space where the generator is uncertain.
That is also the reason why many problems are notsolved (dots in the figure), but even in this challenging
setting our model is more robust.

5

Conclusions

We proposed a simple, fast and robust algorithm to compute shortest paths on Riemannian manifolds learned
from data. Here, standard solvers often fail due to illconditioned Jacobians of the associated BVP. Instead,our method applies a Jacobian-free fixed-point iterative
scheme. The assumption is that the true path can beapproximated smoothly by the predictive GP posterior. This solver makes the Riemannian methods morefeasible since robustly in reasonable time, complex statistical models (Arvanitidis et al., 2016) can be fitted,
as well as distances can be computed in challenging
deep metric scenarios (Arvanitidis et al., 2018).
This has been achieved by analyzing and extending
an existing probabilistic numerical solver (Hennig andHauberg, 2014), turning it from a proof-of-concept
into a principled algorithm. The presented methodthus contributes both to Riemannian methods and toprobabilistic numerics. Further improvements might
be achieved with more complex fixed-point iterations(Ishikawa, 1974), advanced line searches (Mahsereciand Hennig, 2015), adaptive mesh selection (Mazzia
and Trigiante, 2004), and improved model selection(Vaart and Zanten, 2011).

Georgios Arvanitidis† , Søren Hauberg† , Philipp Hennig‡ , Michael Schober?

Acknowledgments
GA is supported by the Danish Center for Big Data
Analytics Driven Innovation. SH was supported bya research grant (15334) from VILLUM FONDEN.
This project has received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (grant agreement no 757360). We gratefully
acknowledge the support of the NVIDIA Corporation
with the donation of the used Titan Xp GPU. PH
gratefully acknowledges financial support by the ERC
StC Action 757275 PANAMA and grant BMBF 01 IS
18 052-B of the German Federal Ministry for Education
and Research.
References
G. Arvanitidis, L. K. Hansen, and S. Hauberg. A locally
adaptive normal distribution. In Advances in Neural
Information Processing Systems (NIPS), 2016.
G. Arvanitidis, L. K. Hansen, and S. Hauberg. Latent space
oddity: on the curvature of deep generative models. In
International Conference on Learning Representations(ICLR), 2018.
U. M. Ascher, R. M. Mattheij, and R. D. Russell. Numerical solution of boundary value problems for ordinary
differential equations, volume 13. Siam, 1994.
N. Bello, A. J. Alkali, and A. Roko. A fixed point iterative method for the solution of two-point boundary value problems for a second order differential equations. Alexandria Engineering Journal, 2017. doi:
10.1016/j.aej.2017.09.010.
F.-X. Briol, C. Oates, M. Girolami, and M. A. Osborne.
Frank-wolfe bayesian quadrature: Probabilistic integration with theoretical guarantees. In Advances in Neural
Information Processing Systems, pages 1162–1170, 2015.
O. A. Chkrebtii, D. A. Campbell, B. Calderhead, and M. A.
Girolami. Bayesian solution uncertainty quantification fordifferential equations. Bayesian Anal., 11(4):1239–1267,
12 2016. doi: 10.1214/16-BA1017.
J. Cockayne, C. Oates, T. Sullivan, and M. Girolami. Probabilistic meshless methods for partial differential equations and bayesian inverse problems. arXiv preprint
arXiv:1605.07811, 2016.
P. Deuflhard. Newton methods for nonlinear problems:
affine invariance and adaptive algorithms. Springer Science & Business Media, 2011.
M. do Carmo. Riemannian Geometry.
(Boston, Mass.). Birkhäuser, 1992.

Mathematics

A. Feragen, F. Lauze, and S. Hauberg. Geodesic exponential kernels: When curvature and linearity conflict. In
Proceedings of the IEEE Conference on Computer Vision
and Pattern Recognition, pages 3032–3042, 2015.
S. Hauberg, O. Freifeld, and M. Black. A Geometric Takeon Metric Learning. In P. Bartlett, F. Pereira, C. Burges,
L. Bottou, and K. Weinberger, editors, Advances in Neural Information Processing Systems (NIPS), pages 2033–
2041. MIT Press, 2012.

S. Hauberg, M. Schober, M. Liptrot, P. Hennig, and A. Feragen. A random riemannian metric for probabilistic
shortest-path tractography. In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015,
volume 18, Munich, Germany, Sept. 2015. Springer.
J.-H. He. Variational iteration method for autonomous
ordinary differential systems. Applied Mathematics andComputation, 114(2):115 – 123, 2000. ISSN 0096-3003.
doi: 10.1016/S0096-3003(99)00104-6.
P. Hennig and S. Hauberg. Probabilistic Solutions to Differential Equations and their Application to Riemannian
Statistics. In Proc. of the 17th int. Conf. on Artificial
Intelligence and Statistics (AISTATS), volume 33. JMLR,
W&CP, 2014.
P. Hennig, M. A. Osborne, and M. Girolami. Probabilistic
numerics and uncertainty in computations. Proceedingsof the Royal Society of London A: Mathematical, Physicaland Engineering Sciences, 471(2179), 2015.
S. Ishikawa. Fixed points by a new iteration method. Proceedings of the American Mathematical Society, 44(1):147–150, 1974.
H. Jafari. A comparison between the variational iteration
method and the successive approximations method. Applied Mathematics Letters, 32:1 – 5, 2014. ISSN 0893-9659.
doi: 10.1016/j.aml.2014.02.004.
G. G. Johnson. Fixed points by mean value iterations.
Proceedings of the American Mathematical Society, 34(1):
193–194, 1972.
M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur. Gaussian Processes and Kernel Methods: A
Review on Connections and Equivalences. ArXiv e-prints,
July 2018.
H. Kersting, T. J. Sullivan, and P. Hennig. Convergence
Rates of Gaussian ODE Filters. ArXiv e-prints, July
2018.
H. P. Kersting and P. Hennig. Active uncertainty calibrationin Bayesian ODE solvers. In Janzing and Ihlers, editors,
Uncertainty in Artificial Intelligence (UAI), volume 32,
2016.
S. Khuri and A. Sayfy. Variational iteration method:
Green’s functions and fixed point iterations perspective.
Applied Mathematics Letters, 32:28 – 34, 2014. ISSN
0893-9659. doi: 10.1016/j.aml.2014.01.006.
G. S. Kimeldorf and G. Wahba. A correspondence between
bayesian estimation on stochastic processes and smoothing by splines. The Annals of Mathematical Statistics, 41(2):495–502, 1970.
D. P. Kingma and M. Welling. Auto-Encoding Variational
Bayes. In Proceedings of the 2nd International Conference
on Learning Representations (ICLR), 2014.
M. Mahsereci and P. Hennig. Probabilistic line searches forstochastic optimization. In Advances in Neural Information Processing Systems, pages 181–189, 2015.
W. R. Mann. Mean value methods in iteration. Proceedings of the American Mathematical Society, 4(3):506–510,
1953.
F. Mazzia and D. Trigiante. A hybrid mesh selection strategy based on conditioning for boundary value ode problems. Numerical Algorithms, 36(2):169–187, 2004.

Fast and Robust Shortest Paths on Manifolds Learned from Data
C. A. Micchelli, Y. Xu, and H. Zhang. Universal kernels.
Journal of Machine Learning Research, 7(Dec):2651–2667,
2006.
C. Oates, J. Cockayne, R. G. Aykroyd, et al. Bayesian
probabilistic numerical methods for industrial processmonitoring. arXiv preprint arXiv:1707.06107, 2017.
E. Picard. Mémoire sur la théorie des équations aux dérivées
partielles et la méthode des approximations successives.
Journal de Mathématiques Pures et Appliqués, 6:145–210,
1890.
C. Rasmussen and C. Williams. Gaussian Processes for
Machine Learning. MIT, 2006.
D. J. Rezende, S. Mohamed, and D. Wierstra. Stochastic
backpropagation and approximate inference in deep generative models. In Proceedings of the 31st International
Conference on Machine Learning, Bejing, China, 2014.
M. Schober, D. Duvenaud, and P. Hennig. Probabilistic
ODE Solvers with Runge-Kutta Means. Advances in
Neural Information Processing Systems (NIPS), 2014.
M. Schober, S. Särkkä, and P. Hennig. A probabilistic model for the numerical solution of initial value
problems. Statistics and Computing, Jan 2018. doi:
10.1007/s11222-017-9798-7.
J. Skilling. Bayesian solution of ordinary differential equations. Maximum Entropy and Bayesian Methods, Seattle,
1991.
O. Teymur, K. Zygalakis, and B. Calderhead. Probabilistic
linear multistep methods. In D. D. Lee, M. Sugiyama,
U. V. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems 29,
pages 4314–4321. Curran Associates, Inc., 2016.
A. Tosi, S. Hauberg, A. Vellido, and N. D. Lawrence. Metrics for Probabilistic Geometries. In The Conference on
Uncertainty in Artificial Intelligence (UAI), July 2014.
A. v. d. Vaart and H. v. Zanten. Information rates of
nonparametric gaussian process methods. Journal of
Machine Learning Research, 12(Jun):2095–2119, 2011.
G. Wahba. Spline models for observational data. Number 59
in CBMS-NSF Regional Conferences series in appliedmathematics. SIAM, 1990.
H. Wendland. Scattered data approximation. Cambridge
University Press, 2004.
X. Xi, F.-X. Briol, and M. Girolami. Bayesian quadrature
for multiple related integrals. In International Conference
on Machine Learning, 2018.