Only Bayes should learn a manifold(on the estimation of differential geometric structure from data)

Søren Hauberg
Section for Cognitive Systems
Technical University of Denmark
sohau@dtu.dk

Abstract
We investigate learning of the differential geometric structure of a data manifold
embedded in a high-dimensional Euclidean space. We first analyze kernel-basedalgorithms and show that under the usual regularizations, non-probabilistic methods
cannot recover the differential geometric structure, but instead find mostly linear
manifolds or spaces equipped with teleports. We repeat the analysis for probabilistic
methods and show that they naturally recover the geometric structure. Fully
exploiting this structure, however, requires the development of stochastic extensions
to classic Riemannian geometry. We take early steps in that regard. Finally, we
partly extend the analysis to models based on neural networks, thereby highlighting
geometric and probabilistic shortcomings of current deep generative models.

Comments on this document are gratefully accepted at sohau@dtu.dk. This is the 2nd revision.

1

Motivation and background

Manifold learning aim to learn a low-dimensional representation of data that reflect the intrinsic
structure of data. Spectral methods seek a low-dimensional embedding of high-dimensional datathat preserve certain aspects of the data. This includes methods such as Isomap [35], Locally linear
embeddings [29], Laplacian eigenmaps [2] and more [32, 9]. Probabilistic methods often view thedata manifold as governed by a latent variable along with a generative model that describe how
the latent manifold is to be embedded in the data space. The common theme is the quest for a
low-dimensional representation that faithfully capture the data.
Ideally, we want an operational representation, i.e. we want to be able to make mathematicallymeaningful calculations with respect to the learned representation. For quantitative data analysis in
the learned representation, a reasonable set of supported “operations” at least include:• Interpolation: given two points, a natural unique interpolating curve that follow the manifold should exist.
• Distances: the distance between two points should be well-defined and reflect the amountof energy required to transform one point to another.
• Measure: the representation should be equipped with a measure under which integration iswell-defined for all points on the manifold.
Depending on which analysis is to be performed in the new representation, one may focus on different
operations. Even if the above operations should be considered elementary, most manifold learning
schemes do not support any of these.
Embedding methods seek a low-dimensional embedding z1:N = {z1 , . . . , zN } of the data x1:N .
These methods fundamentally only describe the data manifold at the points where data is observed
Preprint. Work in progress.

Figure 1:
Reparametrizations illustrated. The left panel shows a“swirling” transformation of Z withthe property that a Gaussian variablewith zero mean and unit covariance,
will have the same distribution after
a reparametrization. The right panel
shows pair-wise distances between
points before and after reparametrization; evidently the geometry of Z is
sensitive to reparametrizations.

and nowhere else. As such, the low-dimensional embedding space is only well-defined at z1:N . It is
common to treat the low-dimensional embedding space as being Euclidean, but this is generally apost hoc assumption with limited grounding in the embedding method. Fundamentally, the learned
representation space is a discrete space that does not lend itself to continuous interpolations. Likewise,
the most natural measure will only assign mass to the points z1:N , and any associated distribution
will be discrete. This is too limited to be considered an operational representation.
Generative models estimate a set of low-dimensional latent variables z1:N along with a suitable
mapping f : Z → X such that f (z) ≈ x. It is, again, common to treat the latent space Z as
being Euclidean. However, this assumption easily lead to arbitrariness. As an example, consider thevariational autoencoder (VAE) [17, 27], which seek a representation in which z1:N follow a unit
Gaussian distribution. Now consider the transformationg(z) = Rθ z,(1.1)

where Rθ is a linear transformation that rotate points by θ(z) = sin(πkzk). This is a smooth
invertible transformation with the property thatz ∼ N (0, I)⇒

g(z) ∼ N (0, I).

(1.2)

Figure 1 illustrate this transformation. If the latent variables z1:N and the generator f is an optimal
VAE, then g(z1:N ) and f ◦ g −1 is equally optimal. Yet, the latent spaces Z and g(Z) are quite
different; Fig. 1 shows the Euclidean distances between zn and g(zn ) for samples drawn from aunit Gaussian. Clearly, the transformed latent space is significantly different from the original space.
As the VAE provides no guarantees as to which latent space is recovered, we must be careful whenrelying on the Euclidean latent space: distances between points are effectively arbitrary and as are
straight-line interpolations. Any analysis relying on vector operations in the latent space are, thus,
arbitrary and positive result should be viewed either as pure luck with little mathematical grounding,or due to unspecified aspects of the model. Ideally, we want a representation space that is invariant tosuch transformations, but current models do not deliver.
In this paper, we consider models where the representation space Z is learned jointly with a smooth
mapping f : Z → X , such that Z is naturally endowed with a Riemannian metric (Sec. 2). This
gives well-defined interpolants, distances and a natural measure. We contribute a detailed analysis of
the case where f is estimated by a kernel method (Sec. 3), and show that even in the case of infinite
noise-free data a non-probabilistic estimate of f cannot recover the true Riemannian structure of
Z. In contrast, we show that probabilistic estimates of f can recover the true Riemannian structure(Sec. 3.2). Fully exploiting this structure, however, require the development of Bayesian extensions
to classic differential geometry (Sec. 4); we contribute elementary results in that regard, but manyquestions remain open. Finally, we partly extend our analysis to the case where f is a neural networkand demonstrate that current deep generative models are lacking elementary properties before they
can learn the Riemannian structure of data manifolds (Sec. 5). Our key finding is that uncertainty
quantification is a prerequisite for learning an operational representation as the usual smoothness
regularization introduce a harmful bias.
Notation. We let Z denote the d-dimensional representation or latent space, which is learned fromdata in the observation space X ≡ RD . Topologically, the latent space Z is assumed to be Euclidean.
Latent points are denoted zn ∈ Z, while corresponding observations are xn ∈ X . The mappingf : Z → X embeds Z in X ; we denote M = f (Z) and assume that M is a Riemannian manifold.
2

2

Riemannian manifolds

A d-dimensional manifold M embedded in RD (d ≤ D) is a topological space in which there exist
a neighborhood around each point x ∈ M that is homeomorphic to Rd [11]. Informally, M is a(usually nonlinear) surface in RD that is locally Euclidean, i.e. it does not self-intersect or otherwise
locally change dimensionality, etc. We assume that we have a d-dimensional parametrization Z ofthe manifold along with a mapping f : Z → X , such that M = f (Z).
P
We first define the inner product between points in RD as hx, x0 i = 1/D i xi x0i . The division by
D ensures that the induced norm remain finite in the limit D → ∞. Now, let z be a d-dimensional
latent point and let ∆1 and ∆2 be infinitesimals, then we can compute their inner product around z in
the data space using Taylor’s Theorem,
hf (z+∆1 )−f (z), f (z+∆2 )−f (z)i
(2.1)
= hf (z)+Jz ∆1 −f (z), f (z)+Jz ∆2 − f (z)i
(2.2)

>
= hJz ∆1 , Jz ∆2 i = 1/D · ∆>(2.3)
1 J z Jz ∆ 2 ,where Jz = ∂z f ∈ RD×d is the Jacobian of f at z. The d × d symmetric positive definite matrix
1/D · J>J , thus defines a local inner product. We denote this matrixz z
Mz = 1/D · J>
(2.4)
z Jz ,
and refer to it as the (pull-back) metric of the manifold. Note that this local inner product is invariant
to reparametrizations of the manifold as it merely correspond to the inner product of X measured
locally on the manifold. Hence it avoids the parametrization issue discussed in the opening section.
Distances & interpolants. The length of a smooth curve in latent space c : [a, b] → Z under the
local inner product is
Z bq
(2.5)
ċ>
L(c) =t Mct ċt dt,awhere ċt = ∂t c(t) is the derivative of the curve. Natural interpolants (geodesics) can then be
defined as length minimizing curves connecting two points. The length of such a curve is a natural
distance measure along the manifold. Unfortunately, minimizing curve length gives rise to a poorly
determined optimization problem as the length of a curve is invariant to its parametrization. The
following proposition provides remedy [11]:Proposition 1. Let c : [a, b] → Z be a smooth curve that (locally) minimizes the “curve energy”
Z
1 b >E(c) =ċ Mct ċt dt,
(2.6)
2 a t
then c has constant velocity and is length-minimizing.
This energy functional is locally uniformly convex and therefore its solution is locally unique.
Standard calculus of variation shows that curves of minimal energy satisfy the following system of
second order differential equations,h1c̈t = − M−1
2(I ⊗ ċ>t )∂ct vec[Mct ] ċt
2 ct
i
>

(2.7)

− ∂ct vec[Mct ] (ċt ⊗ ċt ) ,where vec[·] stacks the columns of a matrix into a vector and ⊗ is the Kronecker product. Such systems can be solved numerically using standard
techniques. Figure 2 gives an example geodesic.
Integration. Given a function h : X → R we can integrate it over a partof the manifold f (Ω), Ω ⊆ Z as [24]ZZph(x)dx =h(f (z)) det(Mz )dz.
(2.8)
f (Ω)

Figure 2: Geodesic interpolation along the manifold (blue) versus along a
straight line (yellow).

Ωp

The quantity det(M) is known as the Riemannian volume measure and is akin to the Jacobiandeterminant in the change of variables theorem. As before, this integration is invariant toreparametrizations of Z as it is performed with respect to the measure of X .
3

Figure 3: (a) Latent points zn ∈ Z for our guiding example. (b) Geodesics for Gaussian kernel ridge
regression. These are pushed away from the data via “teleports”. (c) Geodesics for kernel ridge
regression with a Gaussian+linear kernel. The linear extrapolation implies (almost) linear geodesics.
(d) Geodesics for Gaussian process regression. The uncertainty make geodesics move along the
manifold.

3

Manifold learning with kernels

We now consider data x1:N distributed on a compact d-dimensional Riemannian submanifold M ⊂RD embedded in the data space. We consider a known set of d-dimensional representations z1:Nand estimate the mapping f : Z → X using kernel methods. Note that this manifold is only locallydiffeomorphic to d-dimensional Euclidean space, and it may globally self-intersect1 . For the sake ofanalysis, we assume noise-free data and consider the limit N → ∞. This setting is sufficient to prove
our main point, but the analysis also hold under noise.
Our key question is if we can recover the true Riemannian structure of Z. Methods that fail at thisgiven infinite noise-free data should be avoided. To give an answer, we will study the metric inregions that are near the training data, and in regions that are far away. We formalize this as follows.
Definition 1. For a point z and a dataset Z, the distance between them is dist(z, Z) = inf z̃∈Z kz−z̃k.
Note that this infimum always exist as the element-wise distance is bounded from below by 0.
Definition 2. For a function x = h(z), we define the limitsaway

x −−−→ x̂

near

and x −−−→ x̂

l→∞

(3.1)
l→∞

if for any sequences ẑl such that dist(ẑl , Z) −−−→ ∞, or dist(ẑl , Z) −−−→ 0, respectively, we havel→∞h(ẑl ) −−−→ x̂. Note that these limits are not always defined.
A guiding example. To illustrate our main point, we draw data uniformly on a unit circle andnonlinearly embed it in X = R1000 with added Gaussian noise. We project this data into Z = R2
while keeping the circular structure of data and learn a mapping f from Z to X . Details of thisprocess are in Appendix A. Finally, we compute shortest paths under the pull-back metric; if the true
metric is recovered we should see shortest paths corresponding to circular arcs in Z. Figure 3a show
the latent points in Z.
3.1

The deterministic setting

We now consider learning the mapping f : Z → X using kernel ridge regression [31], i.e.
−1X,fkrr (z∗ ) = k∗,z kz,z + σ 2 I(3.2)

where k is a suitable kernel function, X ∈ RN ×D is the data matrix, and we have used the short-handnotations k∗,z = k(z∗ , z1:N ) ∈ R1×N and kz,z = k(z1:N , z1:N ) ∈ RN ×N . Since we consider
noise-free data, we have σ 2 = 0. The pull-back metric associated with this regression function is−1−1>
Mkrr (z∗ ) = ∂z∗ k∗,z kz,z
XX> kz,z∂z∗ k∗,z.

(3.3)

Assuming a universal kernel [31], then fkrr will correspond to the true mapping where we have datanearwhen N → ∞. Consequently, we recover the true metric where we have data, i.e. Mkrr −−−→ Mtrue .
1

Technically, this render the manifold immersed rather than embedded; the distinction is not important forour purposes.

4Teleports? The behavior away from data depend on the kernel. We first consider the Gaussian
kernel αkRBF (z, z0 ) = θRBF · exp − kz − z0 k2 ,(3.4)
2
but note that similar observations hold for most common stationary kernels. From this, we see thataway

fRBF −−−→ 0

and

away

MRBF −−−→ 0.

(3.5)

To illustrate the geometric implication of this observation, we consider our guiding example. We
compute shortest paths under the pull-back metric; if the true metric is recovered these should becircular arcs in Z. Figure 3b show the recovered geodesics; we see that they systematically shy awayfrom the data and generally do not resemble circular arcs. The explanation is simply that in termsof length-minimization, it is “free” to move through regions where the metric is zero. The result in
Eq. 3.5, thus, implies that geodesics are encouraged to move away from the data. Intuitively, we canthink of regions in Z without data as “teleports” that points can move freely between. This also holdtrue, when the manifold is densely sampled, and consequently geodesics will generally not movealong the data manifold: the manifold geometry is not recovered.
Flat manifolds? These teleports appear because the chosen kernel cause f to extrapolate to a
constant. We now consider a kernel that extrapolate linearly,
kRBF+lin (z, z0 ) = kRBF (z, z0 ) + θlin z> z0 .

(3.6)

Similarly to before, we see that
away

−1
>
fRBF+lin −−−→ θlin z>∗ Zkz,z X = z∗ B,

(3.7)

−1
where B = θlin Zkz,zX ∈ Rd×D . This amounts to linear extrapolation, as expected. When we moveaway from the data, the metric then becomesaway

MRBF+lin −−−→ BB> .

(3.8)

This is a (scaled) Euclidean metric, implying that the learned manifold is flat in regions where we
do not have data. As the pull-back metric measure distances in X , where straight lines are shortest
curves, then geodesics on the learned manifold will be encouraged to go through the flat regions
where data is missing. This is also evident in Fig. 3c that shows the results of our guiding example.
Here we see that geodesics are almost straight lines, implying that the learned manifold did notrecover the structure of the data.
A regularization perspective. The phenomenon can also be understood by considering regressionfunctions f that minimize a local (log) likelihood function [21, 36]L=N
X

wn Ln + λφ[f ],

(3.9)

n=1

where Ln is the loss associated with the nth observation, wn is a weight that decays with the
distance to the point where f is evaluated, and φ is the regularizer. Girosi et al. [13] show that mostregularizers are low-pass filters (thereby avoiding “wiggly” regression functions). Away from the data,
the regularizer dominate the above cost function, implying a strongly low-pass filtered regressionfunction. By Parseval’s theorem [34] this reduce the energy of curves that move away from the data,
which bias geodesics to move away from the data.
An unstable solution? We have seen that the traditional constant and linear extrapolation schemes
imply that we cannot learn the correct geometry: either we introduce teleports or we learn mostly
flat manifolds. We now ask: how should we extrapolate in order to learn the correct geometry? For
geodesics to stay on the manifold, inner products must take large values away from the data.
away

M −−−→ [sufficiently large].
5(3.10)

That is, to ensure that geodesics always stay on the manifold, length-minimization must be penalized
sufficiently for leaving the manifold. We, currently, do not have a tight bound on how large the metricmust be to ensure this, though a loose bound is provided by the radius of the manifold. That is, let
r = sup dist(z, z0 )(3.11)

z,z0 ∈M

denote the largest distance between points on the manifold, then geodesics stay on the manifold ifaway

λmin (Msufficient ) −−−→ r2 .

(3.12)

Here λmin denote the function returning the smallest eigenvalue of a matrix. In differential geometry,it is common to call any locally length-minimizing curve a geodesic. Here we mean the shortest
geodesic. To ensure that all locally length-minimizing curves stay on the manifold, we have no tighter
bound thanaway

λmin (Mideal ) −−−→ ∞.

(3.13)

If the metric must extrapolate to a large matrix, then the Jacobian J must also extrapolate to a large
matrix (since M ∝ J>J). The results of Girosi et al. [13] dictate that regularizing towards functions
with large derivatives imply that the solution is no longer low-pass filtered (i.e. it will “wiggle”).
Regularizing towards large derivatives, thus, go against common wisdom as regression functions withlarge derivatives in regions of little data, will generally not be stable.
Summarizing discussion. Most common choices for estimating f will ensure that the geometry of
M is well-estimated near the data, so the key factor to determine if we can well-estimate the manifold
geometry is the behavior of f away from the data. In these regions, we depend on prior assumptions
on f to determine the geometry. The most common priors are related to the smoothness of f , where
we use the terminology of Girosi et al. [13], i.e. a function is considered “smooth” if its spectrum
is dominated by low-frequency components. We have seen that the smoother assumptions we arewilling to make, the more geodesics are drawn away from the data (“off the manifold”, so to say).
As a specific example, consider a local likelihood (3.9) with regularizer
"
#
2


∂f
φ[f ] = E
= E tr J>J = DE [tr M] .
(3.14)∂z
Bishop [3] has shown that this is (approximately) the regularizer impliedby training under additive noise.2 We see that this common regularizer
imply a “small” metric, when moving away from the data, which in turn
imply short geodesic segments away from the data. Minimization of Figure 4: Smooth intercurve energy (2.6) — or equivalently curve length (2.5) — will, thus, be polation of data imply
biased towards regions of no data as this is where the metric is minimal. that shortest paths measured in data space cut
The general phenomenon is easily understood by considering a data across holes.
manifold with a hole. If the applied regression function is very smooth,then the hole will be interpolated almost linearly, which imply that shortest paths along the manifold
will cross over the hole. In the end, we are, thus, left with a simple choice: either give up on learning
the manifold geometry correctly (by assuming f is very smooth) or give up on stable learning (byassuming f is not very smooth). Neither choice is desirable.
3.2

The Bayesian setting

As a natural probabilistic extension of the previous sections we now let f consist of component-wise
conditionally independent Gaussian processes (GPs) [26],fi (z) ∼ GP(mi (z), k(z, z0 )),∀i = 1, . . . , D.

(3.15)

This is a Gaussian Process latent variable model (GP-LVM) [20]. Here mi and k are the mean andcovariance functions of the ith GP. Note that the posterior mean function coincide with the previously
considered kernel ridge regression model. Like Lawrence [20], we assume the same covariance2

As an example, the mean decoder of a VAE is trained under this regularizer.

6

function across all dimensions to simplify calculations. The geometry of this model was first studied
by Tosi et al. [37].
The pull-back metric M = 1/DJ>J is now a stochastic Riemannian metric since f is stochastic. As
QD
Gaussian variables are closed under differentiation, then J is Gaussian, J ∼ j=1 N (µ(j, :), Σ),and M follows a non-central Wishart distribution [37, 22]
D · M ∼ Wd (D, Σ, Σ−1 E[J]> E[J]).
(3.16)The entire metric by definition follows a generalized Wishart process [38, 16]. A sample pathfrom this process is smooth when the covariance function k is also smooth, and we have a proper
distribution over Riemannian metrics. However, a sample from f gives a manifold that is only locallydiffeomorphic to d-dimensional Euclidean space, and it may globally self-intersect.
Since the metric is a stochastic variable, we cannot apply standard Riemannian geometry to understandthe space Z (curvature is stochastic, geodesics are stochastic, etc). We can, however, inspect theleading moments of the metric1
1
E[M] = E[J>J] = E[J]> E[J] + Σ(3.17)
DD
Σ2ij + Σii Σjj µ>µ> Σµij Σµj+var [Mij ] =+ i 2 .
(3.18)
2
D
D
D
We see that E[M] remain strictly positive, while var [Mij ] = O (1/D) vanishes in the limit D → ∞.
This can equivalently be seen from the central limit theorem. In high dimensions, the metric, thus,
becomes deterministic even if the underlying manifold is stochastic. This is useful as it implies that
we can well-approximate the stochastic metric with a well-understood deterministic metric.
To see if this approach can learn the geometric structure of the data manifold, we again consider theGaussian kernel (3.4). Straight-forward calculations show that
away

near

Σ −−−→ 0andΣ −−−→ αθRBF I,(3.19)
where α and θRBF are the kernel parameters. From this we see that near the data, the expected metric(3.17) coincides with the true pull-back metric of the manifold (as in the deterministic setting),
near 1
E[M] −−−→ E[J]> E[J].
(3.20)
D
In regions of Z where there is no data, we haveaway 1(3.21)
E[M] −−−→ E[J]> E[J] + αθRBF I.
DIf αθRBF is sufficiently large then geodesics will not go through regions of Z where we do nothave data. When data is sampled densely on the manifold, we often estimate large values of α(corresponding to a less smooth manifold), and a large penalty will be payed for “falling off the
manifold”. To validate these observations, we return to our guiding example; Fig. 3d shows thatgeodesics under the expected metric (3.17) actually follow approximately circular arcs. This aligns
with the theoretical analysis and demonstrates that, unlike a deterministic method, a probabilistic
method can actually learn the differential geometric structure of a data manifold. This is, however, noguarantee: we can only accurately learn the manifold geometry when data is sampled sufficiently
dense on the manifold, but now there is hope, whereas deterministic approaches are bound to fail.
Discussion. Deterministic methods can capture local geometry of the data manifold near theobserved data, but they fail to capture the geometry where data is missing. This is not surprising, as
we can generally only learn when we have data. What is, perhaps, more surprising is that if we canestimate the uncertainty of the manifold, then that translate directly into geometric information: if
there is a hole in the manifold, then we can only see it through a lens of uncertainty. Not quantifying
the uncertainty prevents us from seeing holes and boundaries of a data manifold. Uncertainty, thus,
plays the same role as topology in classic geometry, and this must also be estimated from data.

4

Bayesian geometry

As the expected metric can capture the geometry of the data manifold, we seek a better understandingof stochastic Riemannian metrics. We here consider the case of the GP-LVM where f is a smooth GP.
7

4.1

Detour: Euclidean geometry

Before analyzing stochastic Riemannian metrics, consider a stochastic Euclidean metric. Let u, v ∈Rd denote deterministic vectors, and define their stochastic inner producthu, vi = Au>Av = u A>A v,(4.1)

where A ∈ RD×d is a matrix-valued random variable. Under this inner product, the shortest path
between two points is the straight line, so stochasticity does not change our usual intuitions. The
length of this line is, however, stochastic. Its expectation is found by letting ∆ = A(u − v), thenE k∆k2 = E[∆]> E[∆] + tr (cov [∆]) .

(4.2)

The expected distance under a stochastic Euclidean metric grows with both mean and variance of thebasis A. Hence, expected distances are inherently large when the basis (or equivalently the metric)
has large variance.
Example 1 (a Gaussian basis). Let each row of A be drawn independently from N (0, Σ),such that the metric follows a Wishart distribution [22],M = A>A ∼ Wd (D, Σ).

(4.3)

Then the distance from u to v is Nakagami distributed [19]D
2dist(u, v) ∼ Nakagami, Dσu,v,
2(4.4)

and the expected distance is [15]E [dist(u, v)] =

D+1
2
Γ D
2

Γ



√

2σu,v ∝ σu,v ,

(4.5)

2
where σu,v= (u − v)> Σ(u − v).

(4.6)

We see that the expected distance correspond to Mahalanobis’ distance using the inverse
covariance. Also note that scaling Σ also scales the expected distance: very uncertain metrics
imply large expected distances.

4.2

Geodesics

We now return to the geometry of the GP-LVM, and seek to understand shortest paths. Let c :[a, b] → Z denote a deterministic differentiable curve, and let f (c) denote its stochastic embeddingin X . We stress that c is a deterministic curve in Z, while f (c) is a GP in X . The energy (2.6) off (c) is a random quantity and it is natural to consider its expectation with respect to the random
metric. Since the energy integrand is positive, Tonelli’s Theorem tells us that this expected energy is
1Ē(c) ≡ EM [E(f (c))] = EM
2
Z b
1
=
ċ> E [Mct ] ċt dt.
2 a t

"Z
b

#
ċ>
t Mct ċt dt

(4.7)

a

(4.8)

This implies that the curve c with minimal expected energy over the stochastic manifold, is thegeodesic under the deterministic Riemannian metric E [M]. This is exactly the metric considered inSec. 3.2.
We can understand the curve minimizing expected energy in more explicit terms as follows. Let
ut = E[kċt k] and vt = 1 denote two functions over the interval [a, b]; here we use the short-hand8notation kċt k =pċ>t Mct ċt . Then Cauchy-Scwartz’s inequality tells us that2|hu, vi| ≤ kuk2 · kvk2
!2 Z

(4.9)

b

Z b

E[kċt k]2 dt ·≤E[kċt k]dt

Z bdt

(4.10)

E[kċt k]2 dt.

(4.11)

E[kċt k]dt

(4.12)

a

a

a

Z b
= (b − a)
a

Let
"Z
b

#

Z b

kċt kdt =L̄(c) = Eaadenote the expected length of c then
Z b

E[kċt k]2 dt ≥aL̄2 (c).
b−a(4.13)

Equality is achieved when ut and vt are parallel, that is when E[kċt k] is constant. We can always
reparametrize ct to have constant expected speed and achieve equality. Since var [x] = E[x2 ] − E[x]2 ,we see that
Z b
Z b
Z b
E[kċt k]2 dt =E[kċt k2 ]dt −var [kċt k] dt
a

a

a

Z b
= 2Ē(c) −var [kċt k] dt.

(4.14)a

Assuming that the curve has been parametrized to have constant expected speed, we then get
ZL̄2 (c)1 bĒ(c) =+
var [kċt k] dt.
(4.15)
2(b − a) 2 a
Minimizing expected curve energy, thus, does not always minimize the expected curve length. Rather,
this balances the minimization of expected curve length and the minimization of curve variance.
Implications and interpretation. On a deterministic manifold, minimizing curve energy results ina curve of minimal length (by Proposition 1). When the manifold is stochastic, we see that minimizing
expected curve energy does not imply a minimization of expected curve length. In some sense, thisis disappointing. Yet, it is intriguing that minimizing expected energy corresponds to minimizing
a combination of length and variance. Since the expected energy minimizing curve is the geodesic
under the expected Riemannian metric, this also lend itself to easy computation.
Example 2 (a GP prior manifold). There are cases where expected energy and length arestrongly related. Let each dimension of f be a zero-mean GP with a sufficiently smooth
covariance function. Such processes are common for specifying priors. Let the Jacobian of f
at z beJz ∼DYN (0, Σz ),(4.16)

j=1

such that
D · Mz = J>z Jz ∼ Wd (D, Σz ).
The expected energy of a curve c : [a, b] → Z is thenZZ
1 b  >
D b >Ē(c) =
E ċt Mct ċt dt =ċ Σct ċt dt
2 a
2 a t

9

(4.17)(4.18)

>
since ċ>t Mct ċt ∼ W1 (D, ċt Σct ċt ). Following Example 1 we see that the expected lengthof c isZ bq
L̄(c) ∝ċ>
(4.19)t Σct ċt dt.
a

We can then interpret Σz as a Riemannian metric and note that the expressions for curve
length and energy under this metric correspond to Eqs. 4.19 and 4.18. By Proposition 1, wethen have that minimizing expected curve energy also minimize expected curve length.

4.3

Integration

Let Ω ⊆ Z such that f (Ω) ⊆ M, and let h : f (Ω) → R denote a real-valued integratable function
over the specified part of the manifold. We now seek its integral over the entire domain under arandom Riemannian metric M. That is,
Z
Z
p
h̄ =
h(x)dx =h(f (z)) det Mz dz.
(4.20)f (Ω)Ω

Since the metric is stochastic, then so is h̄. We evaluate its expectation asZp
EM [h̄] = EM
h(f (z)) det Mz dzΩ
Z
hpidet Mz dz.
=
h(f (z))E(4.21)(4.22)Ω

That is, the integration is simply performed under the expected volume measure (as expected).
The variance of the integral can be expressed as
 
 
 2
var h̄ = E h̄2 − E h̄ ,(4.23)

where the last term easily is computed from Eq. 4.22. The missing term is
"Z
2 #
p
 2
E h̄ = EM
h(f (z)) det Mz dz,(4.24)Ω

which generally does not permit a closed-form expression.
Since geodesics under the expected metric are well-behaved, it is tempting
to treat the manifold as having this metric. Fromp an integration point-ofview this implies working with the measure det(E[Mz ]) as suggested
by Arvanitidis et al. [1]. The above analysis, however, indicate that it isperhaps more suitable to use the expected measure (4.22).

Figure 5: The expected
measure of the metricassociated with the GPLVM. Notice that the plot
is practically identical to
Fig. 3d.

Example 3 (a GP prior manifold). There are cases where the choice of measure is of lessimportance. To see this, we return to the Gaussian process prior manifold of Example 2.
As before, the metric at z follows a (scaled) Wishart distribution, and Theorem 3.2.15 ofMuirhead’s book [22] tells us that

hp
i p
Γ D+1
2
d
⇔
E
det(DMz ) = 2 det Σz(4.25)
Γ D−d+12
hp
i p
det(Mz ) ∝ det Σz .
(4.26)Ep√
The measure associated with the expected metric is det E [Mz ] = det Σz and we
conclude that for this prior manifold, the two measures are proportional.

10

Example 4 (GP-LVM). Things are not as simple when considering the posterior GP-LVM
manifold. Here the (scaled) metric at z follows a non-central Wishart distribution>
D · M ∼ Wd (D, Σz , Σ−1z E[J] E[J]).

(4.27)

By Theorem 10.3.7 of Muirhead’s book [22] we get that

d
hp
i
p
2 /2 Γ D+2
4
 det Σz
E det(DMz ) =
d/2
D−d+2
πD Γ
4

>
1
D
· 1 F1 − /2, /2, −1/2Σ−1
E[J]E[J]) ,
z(4.28)

where 1 F1 is the confluent hypergeometric function of the first kind. On the other hand, themeasure associated with the expected metric isqp
det E [Mz ] = det (1/D E[Jz ]> E[Jz ] + Σz )(4.29)
and we see that the two measures appear quite different. To understand this difference inpractice, we show the volume measure of the expected metric (4.29) in the backgroundof Fig. 3d. Similarly, we show the expected volume measure (4.28) in Fig. 5. Somewhat
surprisingly, there is no visual difference between the two different measures. This indicatesthat for the GP-LVM, the more simple volume measure associated with the expected metric
may be a good approximation to the expected volume measure.

5

Deep generative models

So far, we have studied kernel based methods due to their ease. The key observations, however,generally hold true, and we now consider neural networks.
5.1

The deterministic setting

The natural case is to estimate f with a (potentially deep) feed-forward neural network. We call thisan autoencoder as these classic methods are the prime example of such an architecture [30], though
we note that other models such as generative adversarial networks [14] also fall within this category.
When we consider the associated pull-back metric, then the same considerations hold true as in the
kernel-based setting (Sec. 3.1). That is, if we regularize towards a smooth f , then geodesics will
naturally cross through holes in the data manifold. To validate this, Fig. 6a shows that the geodesics
of our guiding example are almost straight lines. As before, the lack of uncertainty prevent us fromlearning the manifold topology and geometry.
5.2

The Bayesian setting

In the neural networks literature, (Gaussian) probabilistic mappings f are commonly represented as[23]
f (z) = µ(z) + diag()σ(z), ∼ N (0, ID ),(5.1)

where µ, σ : Rd → RD are neural networks, diag(·) constructs a diagonal matrix from a vector, andID is the D×D identity matrix. This is a key part of variational autoencoders (VAEs) [17, 27] thatrealize the generative model z ∼ N (0, I), x = f (z). From a geometric perspective it is worth notingthat the noise  does not form a smooth process. As such, sample paths from f are not smooth, andone can question the validity of pull-back metrics of this model. If we disregard any such concerns,
then it is easy to show that [1]> (µ) > (σ) 
D · E[Mz ] = J(µ)
Jz
+ J(σ)
Jz ,
(5.2)
z
z
(µ)

(σ)

where Jz and Jz are the Jacobians of µ and σ, respectively. As before, the variance of the metricalso goes to zero when D → ∞, due to the central limit theorem. This can be taken as a hint that the
expected metric is a reasonable geometric structure for Z.
11

Figure 6: Geodesics for various autoencoders. (a) Here f is a smooth feed-forward network. (b) A
naive VAE where both mean and standard deviation of f are smooth feed-forward networks. (c) A
VAE with decaying precision, i.e. the inverse standard deviation of f is a positive RBF network. (d)Same network as (c) but with background color proportional to the expected measure.
A perhaps more sensible way to arrive at Eq. 5.2 is to view Eq. 5.1 as a random projection of thedeterministic manifold spanned byµ(z)f : Rd → R2D ,f (z) =
.
(5.3)
σ(z)
Stacking µ and σ imply that their Jacobians stack as well, such that Eq. 5.2 is the pull-back metricassociated with Eq. 5.3.
Returning to our guiding example, we let µ and σ be smooth feed-forward neural networks. Figure 6bshows that recovered geodesics are almost straight lines, i.e. the model failed to capture the datageometry. This is because σ is a poor proxy for uncertainty [1]. When σ is a feed-forward neuralnetwork, we assume that we can smoothly interpolate the uncertainty estimates recovered at z1:N .
But smooth interpolation of uncertainty is nonsensical.3 To counter this, Arvanitidis et al. [1] modelσ −1 with a positive RBF network [25], which ensure that uncertainty grows away from the data, andFig. 6c shows that it allows us to recover the geometry of the data manifold.
In Example 4, we saw that for the GP-LVM the measure associated with the expected metric waspractically identical to the expected measure, even if their mathematical expressions are quite different.
This does not appear to be the case for the variational autoencoder. Figure 6c shows the measureassociated with the expected metric, Fig. 6d shows the expected measure (computed using sampling).
Unlike for the GP-LVM, we now see significant differences between the two measures. At this stage,
it is unclear which measure is to be preferred from a practical perspective. We see that both measures
exhibit a somewhat arbitrary behavior, which we take as a hint that the RBF network for inverse
variance does not provide an excellent fit to the data.

6

Quantitative summary

Throughout the previous section we have used a simple example
Methodcorr
Hausdorff
GP-LVM
0.996819
0.858773
to illustrate the fundamental bias in deterministic geometry
KRM
0.843259
5.443746
estimation. For completeness, we here quantify these results,
KRM+ridge 0.892995
3.467452
but emphasize that our main objective is to understand thisAE
0.974833
2.825305
bias, rather than to perform empirical studies. Since we know
VAE
0.975200
2.826954
that ground truth geodesics are circular arcs, we can compare
RBF VAE
0.982504
0.788519
these to estimated geodesics. As a first measure of quality,
we compute the correlation between the estimated geodesic Table 1: Correlation between
lengths and the length of the ground truth geodesics. Table 1 length of true and estimatedshow this correlation for each considered model, and individual geodesics, and Hausdorff distances
correlation plots are found in Appendix B. We observe that the between these curves.
GP-LVM achieve an almost perfect correlation, while deterministic kernel methods fare significantly
worse. Individual correlation plots show that for deterministic models short curves provide a better
3
Consider two low-variance temperature readings at the poles of our planet. If we, from this data, interpolate
the temperature at the equator, then a smooth interpolation of the uncertainty would imply a very certain
prediction.

12

correlation than long curves, which is in line with the Riemannian assumption of a locally Euclidean
model. All autoencoder-based models achieve a strong correlation, which is somewhat surprising
given the almost straight geodesics found by deterministic methods.
As a second quality measure, we report the average Hausdorff distance [28] between estimatedand ground truth geodesics (Table 1). We see that the GP-LVM and the VAE with RBF varianceestimation significantly outperform the other methods. This match the visual observations madethroughout the paper.

7

Previous work

Pull-back metrics have been studied in mathematics at least since the seminal work of Gauss [12],and formed the initial foundation of Riemannian geometry. In machine learning, these metrics have
only been studied in few instances. Tosi et al. [37] was the first to give the latent space of the
GP-LVM [20] a geometric foundation. Bishop et al. [4] used the deterministic volume measure of aGenerative Topographic Map (GTM) [5] as a visualization tool; our analysis imply that incorporating
uncertainty should improve such a tool. Recently, several authors have studied the geometry of deep
generative models [33, 7, 18, 1]. Shao et al. [33] and Chen et al. [7] consider pull-back metrics of
VAEs, but only consider the mean of f ; in our terminology they therefore consider autoencoders
rather than variational autoencoders. Shao et al. note that most geodesics in their model are straightlines and speculate that this is because most data manifold are actually flat. Our analysis shows that
this conclusion is most likely incorrect, and that flatness is an artifact of disregarding uncertainty.
Arvanitidis et al. [1] also considered pull-back metrics of variational autoencoders by taking the
expected metric. Here significant curvature is reported, which coincides with intuition.
As an alternative to estimating a geometry that also reflect topology, one can bias the estimatedgeodesics such that they are attracted to data. Chen et al. [7] initialize geodesics by a density
maximizing curve before optimizing curve energy, and in later work Chen et al. [6] force geodesics
to be polygonal curves that interpolate the data. From our perspective, these approaches work arounda deeper more fundamental problem, as it is not clear if these biased curves correspond to geodesics
under any metric (they most surely do not minimize the energy associated with the deterministic
pull-back metric).

8

Concluding remarks

The driving motivation for introducing pull-back metrics in the latent space of a generative model
is to avoid arbitrariness in parametrizing the latent space. This is an important issue if we are tointerpret the latent variables of a fitted model. We have argued that geometry provides a solution
to the issue, but emphasize that this need not be the only solution. We have demonstrated thatmethods that do not quantify their uncertainty cannot, in a meaningful way, capture the geometryof a data manifold. The key issue is that the usual smoothness assumptions imply that holes in thedata manifold are interpolated so smoothly that geodesics are encouraged to pass through the holes
rather than stay on the manifold. Methods that provide reasonable estimates of the uncertainty of theestimated manifold naturally avoid this issue as the uncertainty directly alters the estimated geometry.
We find that uncertainty quantification in manifold learning ends up playing the role of topology inclassic geometry: uncertainty informs us about holes and boundaries in the manifold and provides
us with a global notion of connectivity. Disregarding uncertainty, thus, imply disregarding the mostfundamental aspects manifold learning.
We have provided an extensive analysis of the geometry of Gaussian latent variable models, and havedeveloped an elementary theory for stochastic Riemannian manifolds. This appear to be the first of its
kind. We have seen that minimizing expected curve energy on a stochastic manifold does not imply
minimization of expected curve length, which forces us to reconsider which measure define the mostnatural interpolants. A more formal treatment of this material alongside approximation bounds whenusing expected metrics have recently appeared [10].
Parts of our analysis can be extended to models based on neural networks. This raises two key
issues for future research: 1) since sample paths from deep generative models are not continuous it isperhaps not a good idea to enforce a geometric analysis. It then becomes interesting to determine if
such smoothness can be introduced in deep generative models without sacrificing the computational13

efficiency of the models. 2) Uncertainty is essential for estimating the geometric structure of adata manifold, but current deep generative models provide rather poor estimators of uncertainty. A
heuristic from Arvanitidis et al. [1] seems to work, but more principled methods would be valuable;
Detlefsen et al. [8] provide some early work, but the question largely remain open.
Finally, we repeat the key point of the paper: without uncertainty quantification, we cannot learn thegeometric structure of a data manifold, and any attempt to do so is bound to fail beyond the most
simple examples.
Acknowledgments
The author is grateful to Vagn Lundsgaard Hansen, Martin Jørgensen, Georgios Arvanitidis, Lars
Kai Hansen, David Eklund and Aasa Feragen for enlightening discussions. SH was supported
by a research grant (15334) from VILLUM FONDEN. This project has received funding fromthe European Research Council (ERC) under the European Union’s Horizon 2020 research andinnovation programme (grant agreement no 757360).

References[1] G. Arvanitidis, L. K. Hansen, and S. Hauberg. Latent space oddity: on the curvature of deep generativemodels. In International Conference on Learning Representations (ICLR), 2018.
[2] M. Belkin and P. Niyogi. Laplacian Eigenmaps for Dimensionality Reduction and Data Representation.
Neural Computation, 15(6):1373–1396, June 2003.
[3] C. M. Bishop. Training with noise is equivalent to tikhonov regularization. Neural computation, 7(1):108–116, 1995.
[4] C. M. Bishop, M. Svensen, and C. K. Williams. Magnification factors for the gtm algorithm. 1997.
[5] C. M. Bishop, M. Svensén, and C. K. Williams. Gtm: The generative topographic mapping. Neuralcomputation, 10(1):215–234, 1998.
[6] N. Chen, F. Ferroni, A. Klushyn, A. Paraschos, J. Bayer, and P. van der Smagt. Fast approximate geodesics
for deep generative models. arXiv preprint arXiv:1812.08284, 2018.
[7] N. Chen, A. Klushyn, R. Kurle, X. Jiang, J. Bayer, and P. Smagt. Metrics for deep generative models. In
International Conference on Artificial Intelligence and Statistics, pages 1540–1550, 2018.
[8] N. S. Detlefsen, M. Jørgensen, and S. Hauberg. Reliable training and estimation of variance networks.
arXiv preprint arXiv:1906.03260, 2019.
[9] D. L. Donoho and C. Grimes. Hessian eigenmaps: Locally linear embedding techniques for highdimensional data. Proceedings of the National Academy of Sciences, 100(10):5591–5596, 2003. ISSN
0027-8424.
[10] D. Eklund and S. Hauberg. Expected path length on random manifolds, 2019.
[11] S. Gallot, D. Hulin, and J. Lafontaine. Riemannian geometry, volume 3. Springer, 1990.
[12] C. F. Gauss. Disquisitiones generales circa superficies curvas. Commentationes Societatis Regiae Scientiarum Gottingesis Recentiores, VI:99–146, 1827.
[13] F. Girosi, M. Jones, and T. Poggio. Regularization theory and neural networks architectures. Neuralcomputation, 7(2):219–269, 1995.
[14] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.
Generative adversarial nets. In Advances in neural information processing systems, pages 2672–2680,
2014.
[15] S. Hauberg. The non-central nakagami distribution. 2018.
[16] S. Hauberg. On the geometry of latent variable models. Oberwolfach reports : OWR, (3), 2018.
[17] D. P. Kingma and M. Welling. Auto-Encoding Variational Bayes. In Proceedings of the 2nd International
Conference on Learning Representations (ICLR), 2014.

14

[18] S. Laine. Feature-based metrics for exploring the latent space of generative models, 2018. URL https://openreview.net/forum?id=BJslDBkwG.
[19] D. Laurenson. Nakagami distribution. Indoor Radio Channel Propagation Modelling by Ray Tracing
Techniques, 1994.
[20] N. D. Lawrence. Probabilistic non-linear principal component analysis with gaussian process latent variable
models. Journal of machine learning research, 6(Nov):1783–1816, 2005.
[21] C. Loader. Local Regression and Likelihood. Springer, New York, 1999.
[22] R. J. Muirhead. Aspects of Multivariate Statistical Theory. John Wiley & Sons, 2005.
[23] D. A. Nix and A. S. Weigend. Estimating the mean and variance of the target probability distribution. In
Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94), volume 1, 1994.
[24] X. Pennec. Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements.
Journal of Mathematical Imaging and Vision, 25(1):127–154, July 2006.
[25] Q. Que and M. Belkin. Back to the future: Radial basis function networks revisited. In Artificial Intelligenceand Statistics (AISTATS), 2016.
[26] C. E. Rasmussen and C. Williams. Gaussian Processes for Machine Learning. University Press Group
Limited, 2006.
[27] D. J. Rezende, S. Mohamed, and D. Wierstra. Stochastic backpropagation and variational inference indeep latent gaussian models. In International Conference on Machine Learning, volume 2, 2014.
[28] R. T. Rockafellar and R. J.-B. Wets. Variational analysis, volume 317. Springer Science & BusinessMedia, 2009.
[29] S. T. Roweis and L. K. Saul. Nonlinear dimensionality reduction by locally linear embedding. science,
290(5500):2323–2326, 2000.
[30] D. E. Rumelhart, G. E. Hinton, and R. J. Williams. Learning internal representations by error propagation.
Technical report, California Univ San Diego La Jolla Inst for Cognitive Science, 1985.
[31] B. Schölkopf and A. J. Smola. Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT Press, 2001.
[32] B. Schölkopf, A. Smola, and K.-R. Müller. Kernel principal component analysis. In Advances in Kernel
Methods - Support Vector Learning, pages 327–352, 1999.
[33] H. Shao, A. Kumar, and P. T. Fletcher. The riemannian geometry of deep generative models. arXiv preprint
arXiv:1711.08014, 2017.
[34] J. Y. Stein. Digital signal processing: A computer science perspective. John Wiley & Sons, Inc., 2000.
[35] J. B. Tenenbaum, V. De Silva, and J. C. Langford. A global geometric framework for nonlinear dimensionality reduction. science, 290(5500):2319–2323, 2000.
[36] R. Tibshirani and T. Hastie. Local likelihood estimation. Journal of the American Statistical Association,
82(398):559–567, 1987.
[37] 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.
[38] A. G. Wilson and Z. Ghahramani. Generalised wishart processes. In Uncertainty in Artificial Intelligence(UAI), 2011.

15

A

Experimental details

Data creation: We sample N = 200 numbers {tn }N
n=1 uniformly over the interval [0, 2π], andembed these first in R3 as"#cos(tn )sin(tn )vn =(A.1)
cos(tn ) sin(tn )The data is then generated as v
xn = n + σ,
0

(A.2)

where 0 is a D − 3 dimensional vector of zeros,  is a D dimensional vector drawn from a unit
Gaussian, and σ = 0.1 captures “off manifold”-noise. Here we embed into a D = 1000 dimensional
vector space.
Since the data fundamentally live on a unit circle (which we then nonlinearly embed in RD ), we fix
the latent variables zn to be the first two dimensions of xn ; this correspond to points on the unit circle
with added Gaussian noise. In order to compare different models, we fix the latent variables to havethe same values in all models.
Kernel methods: We first fit a model using Gaussian process (GP) regression from z1:N to x1:N
with hyperparameters estimated using maximum likelihood. We do not update the latent variables
z1:N as is commonly done for the GP-LVM as this would complicate a direct comparison between
different models. The methods based on kernel ridge regression are simple taken as the GP meanfunction. This ensure that the exact same hyperparameters are used in the GP and kernel ridge
regression experiments. Again, this choice was made to simplify the comparison of different
methods.
Neural network methods: We first train a multilayer perceptron from z1:N to x1:N according to
the usual autoencoding criterion. We use this to form our autoencoder models. We keep this mappingas the mean function µ of the variational autoencoders (VAE), and fit the uncertainty σ according
to the usual VAE criterion. As before, we take this restricted approach to ensure that models are ascomparable as possible.

B

Correlation plots

In the main text we quantitatively compare models in two ways. The key idea is to take advantage ofthe fact that we know that ground truth geodesics should be circular arcs. Let cgt denote a groundtruth geodesic, and let cest denote an estimated geodesic. As a first measure of quality, we computethe length of each curve under the model-specific metricZ b
Lgt = L(cgt ) =k∂t cgt kM dt
(B.1)
a

Z b
Lest = L(cest ) =k∂t cest kM dt.

(B.2)

a

Both integrals are evaluated using standard quadrature. In the main text, we report the correlation
between Lgt and Lest for randomly sampled points. More insight into the empirical behavior of thedifferent models can be found be directly plotting Lgt and Lest against each other, which we do in
Fig. 7. Here we see a clear trend that methods with no or meaningless uncertainty quantification all
exhibit a trend that short curves correlate well with the ground truth, but long curves do not. As the
Riemannian model is that we work with locally Euclidean models, this behavior is not too surprising.
We see that with methods with reasonably sensible uncertainty quantification the two lengths have amore well-behaved correlation. It is worth pointing out that that VAE with RBF precision is not aswell behaved as the GP-LVM, which indicate that the variance estimation leaves something to be
desired.
16

Kernel regression with Gaussian kernel

Kernel regression with Gaussian+linear kernel

GP-LVM

Autoencoder

Variational autoencoder

Variational autoencoder with RBF precision

Figure 7: Correlation between length of ground truth geodesics and estimated geodesics for different
models.
As a second measure of quality we consider the Haussdorf distance between cgt and cest . This isdefined as
distH (cgt , cest ) = max{ sup inf kx − yk, sup inf kx − yk}.
x∈cgt y∈cest

We refer to the main text for results and discussion.

17

y∈cest x∈cgt(B.3)