Published as a conference paper at ICLR 2018

L ATENT S PACE O DDITY: ON THE C URVATURE

OF D EEP G ENERATIVE M ODELS

Georgios Arvanitidis, Lars Kai Hansen, Søren Hauberg

Technical University of Denmark, Section for Cognitive Systems
{gear,lkai,sohau}@dtu.dk

A BSTRACT
arXiv:1710.11379v3 [stat.ML] 13 Dec 2021

Deep generative models provide a systematic way to learn nonlinear data distri-
butions through a set of latent variables and a nonlinear “generator” function that
maps latent points into the input space. The nonlinearity of the generator implies
that the latent space gives a distorted view of the input space. Under mild condi-
tions, we show that this distortion can be characterized by a stochastic Rieman-
nian metric, and we demonstrate that distances and interpolants are significantly
improved under this metric. This in turn improves probability distributions, sam-
pling algorithms and clustering in the latent space. Our geometric analysis further
reveals that current generators provide poor variance estimates and we propose a
new generator architecture with vastly improved variance estimates. Results are
demonstrated on convolutional and fully connected variational autoencoders, but
the formalism easily generalizes to other deep generative models.

1 I NTRODUCTION
Deep generative models (Goodfellow et al., 2014; Kingma & Welling, 2014; Rezende et al., 2014)
model the data distribution of observations x ∈ X through corresponding latent variables z ∈ Z
and a stochastic generator function f : Z → X as
x = f (z). (1)
Using reasonably low-dimensional latent variables and highly flexible generator functions allows
these models to efficiently represent a useful distribution over the underlying data manifold. These
approaches have recently attracted a lot of attention, as deep neural networks are suitable generators
which lead to the impressive performance of current variational autoencoders (VAEs) (Kingma &
Welling, 2014) and generative adversarial networks (GANs) (Goodfellow et al., 2014).
Consider the left panel of Fig. 1, which shows the latent representations of digits 0 and 1 from
MNIST under a VAE. Three latent points are highlighted: one point (A) far away from the class
boundary, and two points (B, C) near the boundary, but on opposite sides. Points B and C near the
boundary seem to be very close to each other, while the third is far away from the others. Intuitively,
we would hope that points from the same class (A and B) are closer to each other than to members of
other classes (C), but this is seemingly not the case. In this paper, we argue this seemed conclusion
is incorrect and only due to a misinterpretation of the latent space — in fact points A and B are
closer to each other than to C in the latent representation. Correcting this misinterpretation not
only improves our understanding of generative models, but also improves interpolations, clusterings,
latent probability distributions, sampling algorithms, interpretability and more.
In general, latent space distances lack physical units (making them difficult to interpret) and are sen-
sitive to specifics of the underlying neural nets. It is therefore more robust to consider infinitesimal
distances along the data manifold in the input space. Let z be a latent point and let ∆z1 and ∆z2 be
infinitesimals, then we can compute the squared distance

2 ∂f 
kf (z + ∆z1 ) − f (z + ∆z2 )k = (∆z1 − ∆z2 )| (J|z Jz ) (∆z1 − ∆z2 ), Jz = , (2)
∂z z=z
using Taylor’s Theorem. This implies that the natural distance function in Z changes locally as
it is governed by the local Jacobian. Mathematically, the latent space should not then be seen

1
Published as a conference paper at ICLR 2018

Figure 1: Left: An example of how latent space distances do not reflect actual data distances. Right:
Shortest paths on the surface spanned by the generator do not correspond to straight lines in the
latent space, as is assumed by the Euclidean metric.

as a linear Euclidean space, but rather as a curved space. The right panel of Fig. 1 provides an
example of the implications of this curvature. The figure shows synthetic data from two classes,
and the corresponding latent representation of the data. The background color of the latent space
corresponds to det(J|z Jz ), which can be seen as a measure of the local distortion of the latent
p
space. We interpolate two points from the same class by walking along the connecting straight
line (red); in the right panel, we show points along this straight line which have been mapped by
the generator to the input space. Since the generator defines a surface in the input space, we can
alternatively seek the shortest curve along this surface that connects the two points; this is perhaps
the most natural choice of interpolant. We show this shortest curve in green. From the center panel
it is evident that the natural interpolant is rather different from the straight line. This is due to the
distortion of the latent space, which is the topic of the present paper.

Outline. In Sec. 2 we briefly present the VAE as a representative instance of generative models.
In Sec. 3 we connect generative models with their underlying geometry, and in Sec. 4 we argue that
a stochastic Riemannian metric is naturally induced in the latent space by the generator. This metric
enables us to compute length-minimizing curves and corresponding distances. This analysis, how-
ever, reveals that the traditional variance approximations in VAEs are rather poor and misleading;
we propose a solution in Sec. 4.1. In Sec. 5 we demonstrate how the resulting view of the latent
space improves latent interpolations, gives rise to more meaningful latent distributions, clusterings
and more. We discuss related work in Sec. 6 and conclude the paper with an outlook in Sec. 7.

2 T HE VARIATIONAL AUTOENCODERS ACTING AS THE G ENERATOR

The variational autoencoder (VAE) proposed by Kingma & Welling (2014) is a simple yet powerful
generative model which consists of two parts: (1) an inference network or recognition network
(encoder) learns the latent representation (codes) of the data in the input space X = RD ; and (2) the
generator (decoder) learns how to reconstruct the data from these latent space codes in Z = Rd .
Formally, a prior distribution is defined for the latent representations p(z) = N (0, Id ), and there
exists a mapping function µθ : Z → X that generates a surface in X . Moreover, we assume
that another function σ θ : Z → RD + captures the error (or uncertainty) between the actual data
observation x ∈ X and its reconstruction as x = µθ (z) + σ θ  , where  ∼ N (0, ID ) and 
is the Hadamard (element-wise) product. Then the likelihood is naturally defined as pθ (x | z) =
N (x | µθ (z), ID σ 2θ (z)). The flexible functions µθ and σ θ are usually deep neural networks with
parameters θ.
However, the corresponding posterior distribution pθ (z | x) is unknown, as the marginal likelihood
p(x) is intractable. Hence, the posterior is approximated using a variational distribution qφ (z | x) =
N (z | µφ (x), Id σ 2φ (x)), where the functions µφ : X → Z and σ φ : X → Rd+ are again deep
neural networks with parameters φ. Since the generator (decoder) is a composition of linear maps
and activation functions, its smoothness is based solely on the chosen activation functions.

2
Published as a conference paper at ICLR 2018

The optimal parameters θ and φ are found by maximizing the evidence lower bound (ELBO) of the
marginal likelihood p(x) as
{θ∗ , φ∗ } = argmax Eqφ (z|x) [log(pθ (x|z))] − KL(qφ (z|x)||p(z)), (3)
θ,φ
where the bound follows from Jensen’s inequality. The optimization is based on variations of gra-
dient descent using the reparametrization trick (Kingma & Welling, 2014; Rezende et al., 2014).
Further improvements have been proposed that provide more flexible posterior approximations
(Rezende & Mohamed, 2015; Kingma et al., 2016) or tighter lower bound (Burda et al., 2016).
In this paper, we consider the standard VAE for simplicity. The optimization problem in Eq. 3 is
difficult since poor reconstructions by µθ can be explained by increasing the corresponding variance
σ 2θ . A common trick, which we also follow, is to optimize µθ while keeping σ 2θ constant, and then
finally optimize for the variance σ 2θ .

3 S URFACES AS THE F OUNDATION OF G ENERATIVE M ODELS

Mathematically, a deterministic generative model x = f (z) can be seen as a surface model (Gauss,
1827) if the generator f is sufficiently smooth. Here, we briefly review the basic concepts on sur-
faces, as they form the mathematical foundation of this work.
Intuitively, a surface is a smoothly-connected set of
points embedded in X . When we want to make compu-
tations on a surface, it is often convenient to parametrize
the surface by a low-dimensional (latent) variable z
along with an appropriate function f : Z → X . We
let d = dim(Z) denote the intrinsic dimensionality of
the surface, while D = dim(X ) is the dimensionality of
the input space. If we consider a smooth (latent) curve Figure 2: The Jacobian J of a nonlinear
R1
γ t : [0, 1] → Z, then it has length 0 kγ̇ t kdt, where function f provides pa local basis in the
dγ
γ̇ t = t/dt denotes the velocity of the curve. In prac- input space, while det(J| J) measures
tice, the low-dimensional parametrization Z often lacks the volume of an infinitesimal region.
a principled meaningful metric, so we measure lengths
in input space by mapping the curve through f ,
Z 1 Z 1 
˙
∂f 
Length[f (γ t )] = f (γ ) dt = J γ̇ dt, J = (4)

γt t γt
t

∂z z=γ

0 0 t

where the last step follows from Taylor’s Theorem. This implies that the length of a curve γ t along
the surface can be computed directly in the latent space using the (locally defined) norm
q q q
kJγ γ̇k = (Jγ γ̇)| (Jγ γ̇) = γ̇ | (J|γ Jγ )γ̇ = γ̇ | Mγ γ̇. (5)
Here, Mγ = J|γ Jγ is a symmetric positive definite matrix, which acts akin to a local Mahalanobis
distance measure. This gives rise to the definition of a Riemannian metric, which represents a
smoothly changing inner product structure.
Definition 1. A Riemannian metric M : Z → Rd×d is a smooth function that assigns a symmetric
positive definite matrix to any point in Z.
It should be clear that if the generator function f is sufficiently smooth, then Mγ in Eq. 5 is a
Riemannian metric.
When defining distances across a given surface, it is meaningful to seek the shortest curve connecting
two points. Then a distance can be defined as the length of this curve. The shortest curve connecting
points z0 and z1 is by (trivial) definition
γ (shortest)
t = argmin Length[f (γ t )], γ 0 = z0 , γ 1 = z1 . (6)
γt
A classic result of differential geometry (do Carmo, 1992) is that solutions to this optimization
problem satisfy the following system of ordinary differential equations (ODEs)
| |
γ̈ t = −0.5 · M−1
    
γ t 2 · (γ̇ t ⊗ Id )∂γ t vec Mγ t γ̇ t − ∂γ t vec Mγ t (γ̇ t ⊗ γ̇ t ) (7)
h i
|
= −0.5 · M−1
 
γ t 2 · ∂γ t Mγ t − ∂γ t vec Mγ t (γ̇ t ⊗ γ̇ t ), (8)

3
Published as a conference paper at ICLR 2018

h i
where ∂γ t Mγ t = ∂γ (1) Mγ t , . . . , ∂γ (d) Mγ t with ∂γ (j) Mγ t ∈ Rd×d the partial derivative of Mγ t
t t t
with respect to the j-th component of the curve, vec[·] stacks the columns of a matrix into a vector
2
so ∂γ t vec Mγ t ∈ Rd ×d and ⊗ is the Kronecker product. Note that Eq. 8 enables us to compute
 
the associated Christoffel symbols. For completeness, we provide the derivation of the ODEs in
Appendix A, as well as the steps for computing the Christoffel symbols. Shortest curves can then be
computed by solving the ODEs numerically; our implementation uses bvp5c from Matlab.

4 T HE G EOMETRY OF S TOCHASTIC G ENERATORS

In the previous section, we considered deterministic generators f to provide relevant background
information. We now extend these results to the stochastic case; in particular we consider
f (z) = µ(z) + σ(z)  , µ : Z → X , σ : Z → RD
+ ,  ∼ N (0, ID ). (9)
This is the generator driving VAEs and related models. For our purposes, we will call µ(·) the mean
function and σ 2 (·) the variance function.
Following the discussion from the previous section, it is natural to consider the Riemannian metric
Mz = J|z Jz in the latent space. Since the generator is now stochastic, this metric also becomes
stochastic, which complicates analysis. The following results, however, simplify matters.
Theorem 1. If the stochastic generator in Eq. 9 has mean and variance functions that are at least
twice differentiable, then the expected metric equals
 |    |  
Mz = Ep() [Mz ] = J(µ) z J(µ)
z + J(σ)
z J(σ)
z , (10)
(µ) (σ)
where Jz and Jz are the Jacobian matrices of µ(·) and σ(·).
Remark 1. By Definition 1, the metric tensor must change smoothly, which implies that the Ja-
cobians must be smooth functions as well. This is easily ensured with activation functions for the
neural networks that are C 2 differentiable, e.g. tanh(·), sigmoid(·), and softplus(·).
Theorem 2 (Due to Tosi et al. (2014)). The variance of the metric under the L2 measure vanishes
when the data dimension goes to infinity, i.e. limD→∞ Var (Mz ) = 0.

Theorem 2 suggests that the (deterministic) expected metric Mz is

a good approximation to the underlying stochastic metric when the
data dimension is large. We make this approximation, which allows
us to apply the theory of deterministic generators.
This expected metric has a particularly appealing form, where the
two terms capture the distortion of the mean and the variance func-
(σ) (σ)
tions respectively. In particular, the variance term (Jz )| (Jz )
will be large in regions of the latent space, where the generator has
large variance. This implies that induced distances will be large in
regions of the latent space where the generator is highly uncertain, Figure 3: Example shortest
such that shortest paths will tend to avoid these regions. These paths paths and distances.
will then tend to follow the data in the latent space, c.f. Fig. 3. It is worth stressing, that no learning
is needed to compute this metric: it only consists of terms that can be derived directly through f .

4.1 E NSURING P ROPER G EOMETRY T HROUGH M EANINGFUL VARIANCE F UNCTIONS

Theorem 1 informs us about how the geometry of the generative model depends on both the mean
and the variance of the generator. Assuming successful training of the generator, we can expect to
have good estimates of the geometry in regions near the data. But what happens in regions further
away from the data? In general, the mean function cannot be expected to give useful extrapolations
to such regions, so it is reasonable to require that the generator has high variance in regions that are
not near the data. In practice, the neural net used to represent the variance function is only trained
in regions where data is available, which implies that variance estimates are extrapolated to regions
with no data. As neural nets tend to extrapolate poorly, practical variance estimates tend to be
arbitrarily poor in regions without data.

4
Published as a conference paper at ICLR 2018

Input space Latent space Standard variance estimate Proposed variance estimate
1 0
3 3 3 9
−1
6
0 0 0 −2 0
3
−3
−3 −3 −3 0
data latent repr. −4
-1
-1 0 1 −3 0 3 −3 0 3 −3 0 3

Figure 4: From left to right: training data in X , latent representations in Z, the standard deviation
PD
log( j=1 σj (z)) for the standard variance network, and the proposed solution.

Figure 4 illustrates this problem. The first two panels show the data and its corresponding latent
representations (here both input and latent dimensions are 2 to ease illustration). The third panel
shows the variance function under a standard architecture, deep multilayer perceptron with softplus
nonlinearity for the output layer. It is evident that variance estimates in regions without data are
not representative of either uncertainty or error of the generative process; sometimes variance is
high, sometimes it is low. From a probabilistic modeling point-of-view, this is disheartening. An
informal survey of publicly available VAE implementations also reveals that it is common to enforce
a constant unit variance everywhere; this is further disheartening.
For our purposes, we need well-behaved variance functions to ensure a well-behaved geometry,
but reasonable variance estimates are of general use. Here, as a general strategy, we propose to
model the inverse variance with a network that extrapolates towards zero. This at least ensures that
variances are large in regions without data. Specifically, we model the precision as β ψ (z) = σ21(z) ,
ψ
where all operations are element-wise. Then, we model this precision with a radial basis function
(RBF) neural network (Que & Belkin, 2016). Formally this is written
 
2
β ψ (z) = Wv(z) + ζ, with vk (z) = exp −λk kz − ck k2 , k = 1, . . . , K, (11)
where ψ are all parameters, W ∈ RD×K >0 are the positive weights of the network (positivity ensures
a positive precision), ck and λk are the centers and the bandwidth of the K radial basis functions,
and ζ → 0 is a vector of positive constants to prevent division by zero. It is easy to see that with
this approach the variance of the generator increases with the distance to the centers. The right-
most panel of Fig. 4 shows an estimated variance function, which indeed has the desired property
that variance is large outside the data support. Further, note the increased variance between the two
clusters, which captures that even interpolating between clusters comes with a level of uncertainty. In
Appendix C we also demonstrate that this simple variance model improves the marginal likelihood
p(x) on held-out data.
Training the variance network amounts to fitting the RBF net-
work. Assuming we have already trained the inference network
(Sec. 2), we can encode the training data, and use k-means to
estimate the RBF centers. Then, an estimate for the bandwidths
of each kernel can be computed as
 −2
1 1 X
λk = a kzj − ck k2  (12)
2 |Ck |
zj ∈Ck

where the hyper-parameter a ∈ R+ controls the curvature of

the Riemannian metric, i.e. how fast it changes based on the
uncertainty. Since the mean function of the generator is already
trained, the weights of the RBF can be found using projected
gradient descent to ensure positive weights.
One visualization of the distortion of the latent space
prelative to
the input space is the geometric volume measure det(Mz ), Figure 5: Comparison of (log)
which captures the volume of an infinitesimal area in the input measures of standard (top) and
space. Figure 5 shows this volume measure for both standard proposed (bottom) variances.
variance functions as well as our proposed RBF model. We see that the proposed model captures
the trend of the data, unlike the standard model.

5
Published as a conference paper at ICLR 2018

5 E MPIRICAL R ESULTS
We demonstrate the usefulness of the geometric view of the latent space with several experiments.
Model and implementation details can be found in Appendix D. In all experiments we first train a
VAE and then use the induced Riemannian metric.

5.1 M EANINGFUL D ISTANCES

First we seek to quantify if the induced Rie- Digits Linear Riemannian

mannian distance in the latent space is more
useful than the usual Euclidean distance. {0, 1, 2} 77.57(±0.87)% 94.28(±1.14)%
For this we perform basic k-means cluster- {3, 4, 7} 77.80(±0.91)% 89.54(±1.61)%
ing under the two metrics. We construct 3 {5, 6, 9} 64.93(±0.96)% 81.13(±2.52)%
sets of MNIST digits, using 1000 random
samples for each digit. We train a VAE for Table 1: The F -measure results for k-means.
each set, and then subdivide each into 10 sub-sets, and performed k-means clustering under both
distances. One example result is shown in Fig. 6. Here it is evident that, since the latent points
roughly follow a unit Gaussian, there is little structure to be discovered by the Euclidean k-means,
and consequently it performs poorly. The Riemannian clustering is remarkably accurate. Summary
statistics across all subsets are provided in Table 1, which shows the established F -measure for clus-
tering accuracy. Again, the Riemannian metric significantly improves clustering. This implies that
the underlying Riemannian distance is more useful than its Euclidean counterpart.

Figure 6: The result of k-means comparing the distance measures. For the decision boundaries we
used 7-NN classification.

5.2 I NTERPOLATIONS

Next, we investigate whether the Riemannian metric gives more meaningful interpolations. First,
we train a VAE for the digits 0 and 1 from MNIST. The upper left panel of Fig. 7 shows the latent
space with the Riemannian measure as background color, together with two interpolations. Images
generated by both Riemannian and Euclidean interpolations are shown in the bottom of Fig. 7. The
Euclidean interpolations seem to have a very abrupt change when transitioning from one class to
another. The Riemannian interpolant gives smoother changes in the generated images. The top-
right panel of the figure shows the auto-correlation of images along the interpolants; again we see
a very abrupt change in the Euclidean interpolant, while the Riemannian is significantly smoother.
We also train a convolutional VAE on frames from a video. Figure 8 shows the corresponding latent
space and some sample interpolations. As before, we see more smooth changes in generated images
when we take the Riemannian metric into account.

5.3 L ATENT P ROBABILITY D ISTRIBUTIONS

We have seen strong indications that the Riemannian metric gives a more meaningful view of the la-
tent space, which may also improve probability distributions in the latent space. A relevant candidate
distribution is the locally adaptive normal distribution (LAND) (Arvanitidis et al., 2016)
 
1 2
LAND(z | µ, Σ) ∝ exp − distΣ (z, µ) , (13)
2

6
Published as a conference paper at ICLR 2018

Figure 7: Left: the latent space with example interpolants. Right: auto-correlations of Riemannian
(top) and Euclidean (bottom) samples along the curves of the left panel. Bottom: decoded images
along Euclidean (top rows) and Riemannian (bottom rows) interpolants.

−3

−3 0 3

Figure 8: Left: the latent space and geodesic interpolants. Right: samples comparing Euclidean (top
row) with Riemannian (bottom row) interpolation. Corresponding videos can be found here.

where distΣ is the Riemannian extension of Mahalanobis distance. We fit a mixture of two LANDs
to the MNIST data from Sec. 5.2 alongside a mixture of Euclidean normal distributions. The first
column of Fig. 9 shows the density functions of the two mixture models. Only the Riemannian model
reveals the underlying clusters. We then sample 40 points from each component of these generative
models1 (center column of the figure). We see that the Riemannian model generates high-quality
samples, whereas the Euclidean model generates several samples in regions where the generator
is not trained and therefore produces blurry images. Finally, the right column of Fig. 9 shows all
pairwise distances between the latent points under both Riemannian and Euclidean distances. Again,
we see that the geometric view clearly reveals the underlying clusters.

5.4 R ANDOM WALK ON THE DATA M ANIFOLD

Finally, we consider random walks over the data manifold, which is a common tool for exploring la-
tent spaces. To avoid the walk drifting outside the data support, practical implementations artificially
1
We do not follow common practice and sort samples by their likelihood, as this hides low-quality samples.

7
Published as a conference paper at ICLR 2018

Figure 9: From left to right: the mixture models, generated samples, and pairwise distances. Top
row corresponds to the Riemannian model and bottom row to the Euclidean model.

6
3
4

0 2

Latent repr.
0
−3 Riem annian
Euclidean
−2
End
−3 0 3 − 15 − 10 −5 0

Figure 10: Left: the measure in the latent space. Right: the random walks.

restrict the walk to stay inside the [−1, 1]d hypercube. Here, we consider unrestricted Brownian mo-
tion under both the Euclidean and Riemannian metric. We perform this random walk in the latent
space of the convolutional VAE from Sec. 5.2. Figure 10 shows example walks, while Fig. 11 shows
generated images (video here). While the Euclidean random walk moves freely, the Riemannian
walk stays within the support of the data. This is explained in the left panel of Fig. 10, which shows
that the variance term in the Riemannian metric creates a “wall” around the data, which the random
walk will only rarely cross. These “walls” also force shortest paths to follow the data.

Figure 11: The comparison of the random walks, at the steps 200, 300, 3000, 4000 and 5000.

8
Published as a conference paper at ICLR 2018

6 R ELATED W ORK
Generative models. This unsupervised learning category attracted a lot of attention, especially,
due to the advances on the deep neural networks. We have considered VAEs (Kingma & Welling,
2014; Rezende et al., 2014), but the ideas extend to similar related models. These include exten-
sions that provide more flexible approximate posteriors (Rezende & Mohamed, 2015; Kingma et al.,
2016). GANs (Goodfellow et al., 2014) also fall in this category, as these models have an explicit
generator. While the inference network is not a necessary component in the GAN model, it has been
shown that incorporating it improves overall performance (Donahue et al., 2017; Dumoulin et al.,
2017). The same thoughts hold for approaches that transform the latent space through a sequence of
bijective functions (Dinh et al., 2017)

Geometry in neural networks. Bengio et al. (2013) discuss the importance of geometry in neural
networks as a tool to understand local generalization. For instance, the Jacobian matrix is a measure
of smoothness for a function that interpolates a surface to the given data. This is exactly the implica-
tion in (Rifai et al., 2011), where the norm of the Jacobian acts as a regularizer for the deterministic
autoencoder. Recently, Kumar et al. (2017) used the Jacobian to inject invariances in a classifier.

Riemannian Geometry. Like the present paper, Tosi et al. (2014) derive a suitable Riemannian
metric in Gaussian process (GP) latent variable models (Lawrence, 2005), but the computational
complexity of GPs causes practical concerns. Unlike works that explicitly learn a Riemannian metric
(Hauberg et al., 2012; Peltonen et al., 2004), our metric is fully derived from the generator and
requires no extra learning once the generator is available.

7 D ISCUSSION AND F URTHER E XTENSIONS

The geometric interpretation of representation learning is that the latent space is a compressed and
flattened version of the data manifold. We show that the actual geometry of the data manifold can
be more complex than it first appears.
Here we have initiated the study of proper geometries for generative models. We showed that the
latent space not only provides a low-dimensional representation of the data manifold, but at the same
time, can reveal the underlying geometrical structure. We proposed a new variance network for the
generator, which provides meaningful uncertainty estimates while regularizing the geometry. The
new detailed understanding of the geometry provides us with more relevant distance measures, as
demonstrated by the fact that a k-means clustering, on these distances, is better aligned with the
ground truth label structure than a clustering based on conventional Euclidean distances. We also
found that the new distance measure produces smoother interpolation, and when training Rieman-
nian “LAND” mixture models based on the new geometry, the components aligned much better with
the ground truth group structure. Finally, inspired by the recent interest in sequence generation by
random walks in latent space, we found that geometrically informed random walks stayed on the
manifold for much longer runs than sequences based on Euclidean random walks.
The presented analysis easily extends to sophisticated generative models, where the latent space will
be potentially endowed with more flexible nonlinear structures. This directly implies particularly
interesting geometrical models. An obvious question is: can the geometry of the latent space play
a role while we learn the generative model? Either way, we believe that this geometric perspective
provides a new way of thinking and further interpreting the generative models, while at the same
time it encourages development of new nonlinear models in the representation space.

ACKNOWLEDGMENTS
LKH is supported by Innovation Fund Denmark / the Danish Center for Big Data Analytics Driven
Innovation. SH was supported by a research grant (15334) from VILLUM FONDEN. This project
has received funding from the European Research Council (ERC) under the European Union’s Hori-
zon 2020 research and innovation programme (grant agreement no 757360). We gratefully acknowl-
edge the support of the NVIDIA Corporation with the donation of the used Titan Xp GPU. We thank
Marcelo Hartmann for initiating the discussion related to the Christoffel symbols.

9
Published as a conference paper at ICLR 2018

R EFERENCES
Georgios Arvanitidis, Lars Kai Hansen, and Søren Hauberg. A Locally Adaptive Normal Distribu-
tion. In Advances in Neural Information Processing Systems (NeurIPS), 2016.
Yoshua Bengio, Aaron Courville, and Pascal Vincent. Representation Learning: A Review and New
Perspectives. IEEE Trans. Pattern Anal. Mach. Intell., 35(8):1798–1828, August 2013.
Yuri Burda, Roger B. Grosse, and Ruslan Salakhutdinov. Importance Weighted Autoencoders. In
International Conference on Learning Representations (ICLR), 2016.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using Real NVP. In
International Conference on Learning Representations (ICLR), 2017.
M.P. do Carmo. Riemannian Geometry. Mathematics (Boston, Mass.). Birkhäuser, 1992.
Jeff Donahue, Philipp Krähenbühl, and Trevor Darrell. Adversarial Feature Learning. In Interna-
tional Conference on Learning Representations (ICLR), 2017.
Vincent Dumoulin, Ishmael Belghazi, Ben Poole, Alex Lamb, Martin Arjovsky, Olivier Mastropi-
etro, and Aaron Courville. Adversarially Learned Inference. In International Conference on
Learning Representations (ICLR), 2017.
Carl Friedrich Gauss. Disquisitiones generales circa superficies curvas. Commentationes Societatis
Regiae Scientiarum Gottingesis Recentiores, VI:99–146, 1827.
Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair,
Aaron Courville, and Yoshua Bengio. Generative Adversarial Nets. In Advances in Neural Infor-
mation Processing Systems (NeurIPS), 2014.
Søren Hauberg, Oren Freifeld, and Michael J. Black. A Geometric Take on Metric Learning. In
Advances in Neural Information Processing Systems (NeurIPS) 25, pp. 2033–2041, 2012.
Diederik P Kingma and Max Welling. Auto-Encoding Variational Bayes. In Proceedings of the 2nd
International Conference on Learning Representations (ICLR), 2014.
Diederik P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling.
Improved Variational Inference with Inverse Autoregressive Flow. In Advances in Neural Infor-
mation Processing Systems (NeurIPS), 2016.
Abhishek Kumar, Prasanna Sattigeri, and Tom Fletcher. Improved Semi-supervised Learning
with Gans using Manifold Invariances. In Advances in Neural Information Processing Systems
(NeurIPS). 2017.
Neil Lawrence. Probabilistic non-linear principal component analysis with Gaussian process latent
variable models. Journal of machine learning research, 6(Nov):1783–1816, 2005.
Jaakko. Peltonen, Arto Klami, and Samuel Kaski. Improved learning of riemannian metrics for
exploratory analysis. Neural Networks, 17(8):1087–1100, 2004.
Qichao Que and Mikhail Belkin. Back to the future: Radial basis function networks revisited. In
Artificial Intelligence and Statistics (AISTATS), 2016.
Danilo Rezende and Shakir Mohamed. Variational Inference with Normalizing Flows. In Proceed-
ings of the 32nd International Conference on Machine Learning (ICML), 2015.
Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic Backpropagation and
Approximate Inference in Deep Generative Models. In Proceedings of the 31st International
Conference on Machine Learning (ICML), 2014.
Salah Rifai, Pascal Vincent, Xavier Muller, Xavier Glorot, and Yoshua Bengio. Contractive Auto-
Encoders: Explicit Invariance During Feature Extraction. In Proceedings of the 28th International
Conference on Machine Learning (ICML), 2011.
Alessandra Tosi, Søren Hauberg, Alfredo Vellido, and Neil D. Lawrence. Metrics for Probabilistic
Geometries. In The Conference on Uncertainty in Artificial Intelligence (UAI), July 2014.

10
Published as a conference paper at ICLR 2018

A T HE D ERIVATION OF THE G EODESIC D IFFERENTIAL E QUATION

The shortest path between two points x, y ∈ M on a Riemannian manifold M is found by optimiz-
ing the functional
Z 1q
(shortest)
γt = argmin hγ̇ t , Mγ t γ̇ t idt, γ(0) = x, γ(1) = y (14)
γt 0

where γ t : [0, 1] → M and γ̇ t = ∂γ ∂t . The minima of this problem can be found instead by
t

optimizing the curve energy (do Carmo, 1992), so the functional becomes
1 1
Z
(shortest)
γt = argmin hγ̇ t , Mγ t γ̇ t idt, γ(0) = x, γ(1) = y. (15)
γt 2 0
The inner product can be written explicitly as
d X
d
X (i) (j) |
· Mγ(ij)

L(γ t , γ̇ t , Mγ t ) = hγ̇ t , Mγ t γ̇ t i = γ̇t · γ̇t t
= vec Mγ t (γ̇ t ⊗ γ̇ t ) (16)
i=1 j=1

where the index in the parenthesis represents the corresponding element in the vector or matrix. In
the derivation the ⊗ is the usual Kronecker product and the vec[·] stacks the column of a matrix into
  2
a vector so vec Mγ t ∈ Rd . We find the minimizers by the Euler-Lagrange equation
∂L ∂ ∂L
= (17)
∂γ t ∂t ∂ γ̇ t
where
∂ ∂hγ̇ t , Mγ t γ̇ t i
 
∂ ∂L ∂
∂Mγ t
= = 2 · Mγ t γ̇ t = 2 γ̇ t + Mγ t γ̈ t . (18)
∂t ∂ γ̇ t ∂t ∂ γ̇ t ∂t ∂t
Since the term
| |
∂Mγ(11) ∂Mγ(1D)
 
∂Mγ(11) ∂Mγ(1D)
 
t
··· t
t
∂γ t γ̇ t ··· t
∂γ t γ̇ t
∂t ∂t |
  (2D) |

(21)
∂Mγ t (2D)
∂Mγ t ∂Mγ(21) ∂Mγ

∂Mγ t  ···   t
γ̇ t ··· γ̇ t

∂t ∂t   ∂γ t ∂γ t

= = , (19)

∂t .. .. ..   .. .. .. 

 . . .   . . . 
∂Mγ(D1) ∂Mγ(DD)
 | | 
∂Mγ(D1) ∂Mγ(DD)
∂t ··· ∂t
t t
γ̇ t ··· t
γ̇ t
∂γ t ∂γ t

we can write the right hand side of the Eq. 17 as

"   #
∂ ∂L | ∂vec Mγ t
= 2 (γ̇ t ⊗ Id ) γ̇ t + Mγ t γ̈ t . (20)
∂t ∂ γ̇ t ∂γ t

The left hand side term of the Eq. 17 is equal to

 |
∂L ∂  |
∂vec Mγ t
= vec Mγ t (γ̇ t ⊗ γ̇ t ) = (γ̇ t ⊗ γ̇ t ). (21)
∂γ t ∂γ t ∂γ t

The final system of 2nd order ordinary differential equations is

"    | #
1 −1 | ∂vec Mγ t ∂vec Mγ t
γ̈ t = − Mγ t 2 · (γ̇ t ⊗ Id ) γ̇ t − (γ̇ t ⊗ γ̇ t ) . (22)
2 ∂γ t ∂γ t
h i
∂vec[Mγ t ]
Note that (γ̇ |t ⊗Id ) ∂γ t γ̇ t = ∂γ t Mγ t (γ̇ t ⊗ γ̇ t ), where ∂γ t Mγ t = ∂γ (1) Mγ t , . . . , ∂γ (d) Mγ t
t t
d×d
with ∂γ (j) Mγ t ∈ R the partial derivative of Mγ t with respect to the j-th component of the curve.
t
γt ∂vec[M ]
This fact and the notation ∂γ t vec[Mγ t ] = ∂γ t allows us to rewrite Eq. 22 as
1 h | i
γ̈ t = − M−1

2 · ∂γt Mγt − ∂γt vec M γt (γ̇ t ⊗ γ̇ t ). (23)
2 γt

11
Published as a conference paper at ICLR 2018

The classical geodesic equation in Riemannian geometry (do Carmo, 1992) is written as
d X
d d
(k) (i) (j)
X X
γ̈t = Γkij γ̇t γ̇t , Γkij = mkl (∂i mjl + ∂j mli − ∂l mij ), (24)
i=1 j=1 l=1

(ij)
where Γkij are the Christoffel symbols with mij = Mγ−1
t
the i, j element of the inverse metric
(ij)
and ∂k mij = ∂γ (k) Mγ t the partial derivative of the metric element i, j with respect to the k-th
t
component of the curve γ t . Note that these Christoffel symbols are symmetric i.e. Γkij = Γkji .
Even if the derived ODE system in Eq. 23 is equivalent to Eq.h 24, the corresponding coefficients are
−1
 | i d×d2
not equal. In particular, let the coefficient matrix C = Mγ t 2·∂γ t Mγ t −∂γ t vec Mγ t ∈R .
We can easily see that C (ij) 6= C (ji) , however, comparing Eq. 23 with Eq. 24 shows that due to
symmetries the following relation holds between our coefficients C and the true Christoffel symbols
C (kp) + C (kq )
Γkij = Γkji = , with p = d(i − 1) + j and q = d(j − 1) + i. (25)
2

B T HE D ERIVATION OF THE R IEMANNIAN M ETRIC

The proof of Theorem 1.

Proof. As we introduced in Eq. 9 the stochastic generator is

f (z) = µ(z) + σ(z)  , µ : Z → X , σ : Z → RD
+ ,  ∼ N (0, ID ). (26)
Thus, we can compute the corresponding Jacobian as follows
∂fz(1) ∂fz(1) ∂fz(1)
 
∂z 1 ∂z 2
· · · ∂z d
 ∂f (2) ∂fz(2) ∂fz(2) 

∂f (z) 
 ∂z1
z
∂z2 ··· ∂zd
= Jz =  (27)

∂z .. .. .. .. 

 . . . .


(D) (D) (D)
∂fz ∂fz ∂fz
∂z1 ∂z2 · · · ∂zd D×d
∂µ(1) (1) (1)
 
∂µ ∂µ
∂z1
z
∂z2
z
··· ∂zd
z
 ∂µ (2) (2) (2)
∂µ ∂µ


 ∂z1
z
∂z2
z
· · · ∂zd
z 
= + [S1 , S2 , · · · , Sd ]D×d , (28)

.. .. .. .. 

. . . .
 | {z }
B
 
∂µ(D) ∂µ (D)
∂µ (D)

∂z1
z
∂z2
z
··· ∂zd
z
D×d
| {z }
A
∂σz(1)
 
0 ··· 0
 ∂zi ∂σz(2)

 0
∂zi ··· 0 
where Si =  . , i = 1, . . . , d (29)
 
 . .. .. .. 
 . . . .


∂σz(D)
0 0 ··· ∂zi D×D

and the resulting “random” metric in the latent space is Mz = J|z Jz . The randomness is due to the
random variable , and thus, we can compute the expectation
Mz = Ep() [Mz ] = Ep() [(A + B)| (A + B)] = Ep() [A| A + A| B + B| A + B| B]. (30)

Using the linearity of expectation we get that

Ep() [A| B] = Ep() [A| [S1 , S2 , · · · , Sd ]] = A| [S1

Ep()
 :.0. . , 0] = 0
[],
 (31)

12
Published as a conference paper at ICLR 2018

because Ep() [] = 0. The other term

 
| S1

 | S2  
Ep() [B| B] = Ep()  . [S , S , · · · , S ] (32)
 
1 2 d

 .. 


| Sd d×D
 |
 S1 S1  | S1 S2  · · · | S1 Sd 

 | S2 S1  | S2 S2  · · · | S2 Sd  
= Ep()   .. .. .. ..  (33)
 . . . . 
| Sd S1  | Sd S2  · · · | Sd Sd 
∂σz(1)
  
1 ∂zj
∂σ (2)
 ! 
(1) (2) (D)
|
 ∂σz ∂σz ∂σz  2 ∂zzj 
with Ep() [ Si Sj ] = Ep()  1 , 2 , · · · , D (34)
  
∂zi ∂zi ∂zi
 .. 
.
  
  
∂σz(D)
D ∂zj
" ! ! !#
(1) (1) (2) (2) (D) (D)
∂σz ∂σz ∂σz ∂σz ∂σz ∂σz
= Ep() 21 + 22 + · · · 2D
∂zi ∂zj ∂zi ∂zj ∂zi ∂zj
(35)
|
= diag(Si ) diag(Sj ), (36)

because Ep() [2i ] = 1, ∀i = 1, . . . , D.

(µ)
The matrix A = Jz and for the variance network
(1) (1)
∂σ (1)
 
∂σ z ∂σ z
∂z1 ∂z2 ··· z
∂zd
∂σ (2) ∂σ (2) ∂σ (2)
 
 z
∂z1
z
∂z2 ··· z
∂zd

J(σ) = (37)
 
z .. .. .. .. 

 . . . .


∂σ (D) ∂σ (D) ∂σ (D)
z
∂z1
z
∂z2 ··· z
∂zd
 |
(σ) (σ)
it is easy to see that Ep() [B| B] = Jz Jz . So the expectation of the induced Riemannian
metric in the latent space by the generator is
 |  |
M̄z = J(µ)
z Jz(µ) + J(σ)
z J(σ)
z (38)

which concludes the proof.

C I NFLUENCE OF VARIANCE ON THE M ARGINAL L IKELIHOOD

We trained a VAE on the digits 0 and 1 of the MNIST scaled to [−1, 1]. We randomly split the data
to 90% training and 10% test data, ensuring balanced classes. First, we only trained the encoder and
the mean function of the decoder. Then, keeping these fixed, we trained two variance functions: one
based on standard deep neural network architecture, and the other using our proposed RBF model.
Clearly, we have two generators with the same mean function, but different variance functions.
Below we present the architectures for the standard neural networks. For the RBF model we used
32 centers and a = 1.

Encoder/Decoder Layer 1 Layer 2 Layer 3

µφ 64, (softplus) 32, (softplus) d, (linear)
σφ 64, (softplus) 32, (softplus) d, (softplus)
µθ 32, (softplus) 64, (softplus) D, (tanh)
σθ 32, (softplus) 64, (softplus) D, (softplus)

13
Published as a conference paper at ICLR 2018

The numbers corresponds to the layer size together with the activation function in parenthesis. Fur-
ther, the mean and the variance functions share the weights of the first layer. The input space
dimension is D = 784. Then, we computed the marginal likelihood p(x) of the test data using
Monte Carlo as:

Z
1X
p(x) = p(x|z)p(z)dz ' p(x|zs ), zs ∼ p(z) (39)
Z S s=1
using S = 10000 samples. The generator with the standard variance function achieved -68.25 mean
log-marginal likelihood, while our proposed model -50.34, where the higher the better.
The reason why the proposed RBF model performs better can be easily analyzed. The marginal
likelihood under the Monte Carlo estimation is, essentially, a large Gaussian mixture model with
equal weights S1 . Each mixture
component is defined by the generator through the likelihood
p(x|z) = N x | µθ (z), ID σ 2θ (z) . Considering the variance term, the standard neural network
approach is trained on the given data points and the corresponding latent codes. Unfortunately,
its behavior is arbitrary in regions where there are not any encoded data. On the other hand our
proposed model assigns large variance to these regions, while on the regions where we have latent
codes its behavior will be approximately the same with the standard neural network. This implies
that the resulting marginal likelihood p(x) for the two models are highly similar in regions of high
data density, but significantly different elsewhere. The RBF variance model ensures that mixture
components in these regions have high variance, whereas the standard architecture assign arbitrary
variance. Consequently, the RBF-based p(x) assigns minimal density to regions with no data, and,
thus, attains higher marginal likelihood elsewhere.

D I MPLEMENTATION D ETAILS FOR THE E XPERIMENTS

Algorithm 1 The training of a VAE that ensures geometry

Output: the estimated parameters of the neural networks θ, φ, ψ
1: Train the µφ , σ φ , µθ as in Kingma & Welling (2014), keeping σ ψ fixed.
2: Train the σ ψ as explained in Sec. 4.1.

Details for Experiments 5.1, 5.2 & 5.3. The pixel values of the images are scaled to the interval
[0, 1]. We use for the functions µφ , σ φ , µθ multilayer perceptron (MLP) deep neural networks, and
for the β ψ the proposed RBF model with 64 centers, so W ∈ RD×64 and the parameter a of Eq. 12
is set to 2. We used L2 regularization with parameter equal to 1e−5 .

Encoder/Decoder Layer 1 Layer 2 Layer 3

µφ 64, (tanh) 32, (tanh) d, (linear)
σφ 64, (tanh) 32, (tanh) d, (softplus)
µθ 32, (tanh) 64, (tanh) D, (sigmoid)

The number corresponds to the size of the layer, and in the parenthesis the activation function.
For the encoder, the mean and the variance functions share the weights of the Layer 1. The input
space dimension D = 784. After the training, the geodesics can be computed by solving Eq. 7
numerically. The LAND mixture model is fitted as explained in (Arvanitidis et al., 2016).

Details for Experiments 5.4. In this experiment we used Convolutional Variational Auto-
Encoders. The pixel values of the images are scaled to the interval [0, 1]. For the β ψ we used
the proposed RBF model with 64 centers and the parameter a of Eq. 12 is set to 2.
Considering the variance network during the decoding stage, the RBF generates an image, which
represents intuitively the total variance of each pixel for the decoded final image, but in an initial
sub-sampled version. Afterwards, this image is passed through a sequence of deconvolution layers,
and at the end will represent the variance of every pixel for each RGB channel. However, it is critical
that the weights of the filters must be clipped during the training to R+ to ensure positive variance.

14
Published as a conference paper at ICLR 2018

Encoder Layer 1 (Conv) Layer 2 (Conv) Layer 3 (MLP) Layer 4 (MLP)

µφ 32, 3, 2, (tanh) 32, 3, 2, (tanh) 1024, (tanh) d, (linear)
σφ 32, 3, 2, (tanh) 32, 3, 2, (tanh) 1024, (tanh) d, (softplus)

For the convolutional and deconvolutional layers, the first number is the number of applied filters,
the second is the kernel size, and third is the stride. Also, for the encoder, the mean and the variance
functions share the convolutional layers. We used L2 regularization with parameter equal to 1e−5 .

Decoder L. 1 (MLP) L. 2 (MLP) L. 3 (DE) L. 4 (DE) L. 5 (DE) L.6 (CO)

µθ 1024, (t) D/4, (t) 32, 3, 2, (t) 32, 3, 2, (t) 3, 3, 1, (t) 3, 3, 1, (s)

For the decoder, the acronyms (DE) = Deconvolution, (CO) = Convolution and (t), (s) stand for
tanh and sigmoid, respectively. Also, D = width × height × channels of the images, in our
case 64,64,3. For all the convolutions and deconvolutions, the padding is set to same. We used L2
regularization with parameter equal to 1e−5 .

Decoder Layer 1 (RBF) Layer 2 (Deconv) Layer 3 (Conv)

βψ W ∈ R(D/2)×64 1, 3, 2 (linear) 3, 3, 1 (linear)

The Brownian motion over the Riemannian manifold in the latent space is presented in Alg. 2.

Algorithm 2 Brownian motion on a Riemannian manifold

Input: the starting point z ∈ Rd×1 , stepsize s, number of steps Ns , the metric tensor M(·).
Output: the random steps Z ∈ RNs ×d .
1: for n = 0 to Ns do
2: L, U = eig (M(z)), (L: eigenvalues, U: eigenvectors)
− 21
3: v = UL ,  ∼ N (0, Id )
4: z=z+s·v
5: Z(n, :) = z
6: end for

15