Maximum likelihood estimation of Riemannian
metrics from Euclidean data
Georgios Arvanitidis, Lars Kai Hansen, and Søren Hauberg
Section for Cognitive Systems, Technical University of Denmark,Richard Petersens Plads, B321, 2800 Kgs. Lyngby, Denmark.
{gear, lkai, sohau}@dtu.dk
Abstract. Euclidean data often exhibit a nonlinear behavior, which
may be modeled by assuming the data is distributed near a nonlinear
submanifold in the data space. One approach to find such a manifold
is to estimate a Riemannian metric that locally models the given data.
Data distributions with respect to this metric will then tend to follow
the nonlinear structure of the data. In practice, the learned metric relyon parameters that are hand-tuned for a given task. We propose to estimate such parameters by maximizing the data likelihood under theassumed distribution. This is complicated by two issues: (1) a changeof parameters imply a change of measure such that different likelihoods
are incomparable; (2) some choice of parameters renders the numericalcalculation of distances and geodesics unstable such that likelihoods cannot be evaluated. As a practical solution, we propose to (1) re-normalize
likelihoods with respect to the usual Lebesgue measure of the data space,and (2) to bound the likelihood when its exact value is unattainable. We
provide practical algorithms for these ideas and illustrate their use on
synthetic data, images of digits and faces, as well as signals extracted
from EEG scalp measurements.
Keywords: manifold learning, metric learning, statistics on manifolds.

1

Introduction

The “manifold assumption” is often applied in machine learning research toexpress that data is believed to lie near a (nonlinear) submanifold embeddedin the data space. Such an assumption finds uses e.g. in dynamical or periodic
systems, and in many problems with a smooth behavior. When the manifold
structure is known a priori it can be incorporated into the problem specification,
but unfortunately such structure is often not known. In these cases it is necessary
to estimate the manifold structure from the observed data, a process known as
manifold learning. In this work, we approach manifold learning geometricallyby estimating a Riemannian metric that captures local behavior of the data,and probabilistically by estimating unknown parameters of the metric using
maximum likelihood. First we set the stage with background information on
manifold learning (Sec. 1.1) and geometry (Sec. 1.2), followed by an expositionof our model (Sec. 2) and the proposed maximum likelihood scheme (Sec. 3).
Finally results are presented (Sec. 4) and discussed (Sec. 5).

1.1

Background and related work

Given observations x1:N = {x1 , . . . , xN } in RD , the key task in manifold learning
is to estimate a data representation that reflect the nonlinear structure of theoriginal data. The intuition behind most methods for this was phrased by Saul
& Roweis [18] as “Think Globally, Fit Locally”, practically meaning that locallylinear models are fitted to all points in data space and these then are merged toa global representation (details depend on the method).
The Isomap method [20] famously replace Euclidean distances with geodesic
distances defined on a neighborhood graph and then embed the data in a lower
dimensional space where Euclidean distances approximate the geodesic counterparts. While this approach is popular, its discrete nature only describes theobserved data points and consequently cannot be used to develop probabilistic
generative models. Similar comments hold for other graph-based methods [18,2].
As a smooth alternative, Lawrence [16] proposed a probabilistic extension ofstandard surface models by assuming that each dimension of the data is described
as xd = fd (z), where z is a low-dimensional latent variable and fd is a Gaussian
process. The latent variables then provide a low-dimensional parametrization that
capture the manifold strucure. Tosi et al. [21] give this a geometric interpretation
by deriving the distribution of the induced Riemannian pull-back metric andshow how geodesics can be computed under this uncertain metric.
Often manifold learning is viewed as a form of dimensionality reduction,
but this need not be the case. Hauberg et al. [12] suggest to model the local
behavior of the data manifold via a locally-defined Riemannian metric, which isconstructed by interpolating a set of pre-trained metric tensors at a few selectpoints in data space. Once a Riemannian metric is available existing tools canbe used for dimensionality reduction [22,11,8], mixture modeling [1,19], tracking[13,14], hypothesis testing [17], transfer learning [9] and more. Our approach
follow this line of work.
1.2

The basics of Riemannian geometry

For completeness we start with an informal review of Riemannian manifolds, butrefer the reader to standard text books [5] for a more detailed exposition.
Definition 1. A smooth manifold M together with a Riemannian metric M :M → RD×D and M  0 is called a Riemannian manifold. The Riemannian
metric M encodes a smoothly changing inner product hu, M(x)vi on the tangent
space u, v ∈ Tx M of each point x ∈ M.
Since the Riemannian metric M(x) acts on tangent vectors it may be interpreted
as a standard Mahalanobis metric restricted to an infinitesimal region around x.
This local inner product is a suitable model for capturing local behavior of data,
i.e. manifold learning. Shortest paths (geodesics) are then length-minimizing
curves connecting two points x, y ∈ M, i.e.
Z 1pγ̂ = argmin
hγ 0 (t), M(γ(t))γ 0 (t)idt, s.t. γ(0) = x, γ(1) = y.
(1)
γ

0

Here M(γ(t)) is the metric tensor at γ(t), and the tangent vector γ 0 denotes thederivative (velocity) of γ. The distance between x and y is defined as the lengthof the geodesic. Geodesic can be found as the solution to a system of 2nd order
ordinary differential equations (ODEs):

|
1 −1
∂vec[M(γ(t))]
00
γ (t) = − M (γ(t))(γ 0 (t) ⊗ γ 0 (t))
(2)
2
∂γ(t)
subject to γ(0) = x, γ(1) = y. Here vec[·] stacks the columns of a matrix into avector and ⊗ is the Kronecker product.
This differential equation allows us to define basic operations on the manifold.
The exponential map at a point x takes a tangent vector v ∈ Tx M to y =Expx (v) ∈ M such that the curve γ(t) = Expx (t · v) is a geodesic originating
at x with initial velocity v and length kvk. The inverse mapping, which takes y
to Tx M is known as the logarithm map and is denoted Logx (y). By definition
kLogx (y)k corresponds to the geodesic distance from x to y. The exponential
and the logarithmic map can be computed by solving Eq. 2 numerically, as an
initial value problem or a boundary value problem respectively.

2

A locally adaptive normal distribution

We have previously provided a simple nonparametric manifold learning schemethat conceptually mimics a local principal component analysis [1]. At eachpoint x ∈ RD a local covariance matrix is computed and its inverse thenspecify a local metric. For computational efficiency and to prevent overfitting
we restrict ourselves to diagonal covariances
!−1
N
X
2
wn (x)(xnd − xd ) + ρ ,(3)
Mdd (x) =n=1
2

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

!
.

(4)

Here the subscript d is the dimension, n corresponds
Fig. 1: Example geodesics.
to the given data, and ρ is a regularization parameterto avoid singular covariances. The weight-function wn (x) changes smoothly suchthat the resulting metric is Riemannian. It is easy to see that if x is outsideof the support of the data, then the metric tensor is large. Thus, geodesics are“pulled” towards the data where the metric is small (see Fig. 1).
The weight-function wn (x) depends on a parameter σ that effectively determine the size of the neighborhood used to define the data manifold. Small
values of σ gives a manifold with high curvature, while a large σ gives an almost
flat manifold. The main contribution of this paper is a systematic approach todetermine this parameter.
For a given metric (and hence σ), we can estimate data distributions with respect to this metric. We consider Riemannian normal distributions [17]

pM (x | µ, Σ) =11
exp − d2Σ (x, µ) ,C2x ∈ M (5)and mixtures thereof. Here M denote the manifold induced by the learned metric, µ and Σ are the mean andcovariance, and d2Σ (x, µ) = hLogµ (x), Σ −1 Logµ (x)i.
The normalization constant C is by definition

Z
Fig. 2: Example of the
1 2
locallyadaptive normalC(µ, Σ) =exp − dΣ (x, µ) dM(x), (6)2M
distribution (LAND).
where dM(x) denotes the measure induced by the Riemannian metric. Note that
this measure depends in σ. Figure 2 show an example of the resulting distribution
under the proposed metric. As the distribution adapts locally to the data we coin
it a locally adaptive normal distribution (LAND).
Assuming that the data are generated from a distribution qM then commonly
the mean µ and covariance Σ are estimated with intrinsic least squares (ILS)
Intrinsic leastsquares (ILS)Zµ̂ = argmin
d2 (µ, x)qM (x)dM(x),
(7)
µ∈M
M
Z
Σ̂ =
Logµ̂ (x)Logµ̂ (x)| pM (x)dM(x), (8)Tµ̂ MMaximum
likelihood (ML)where d2 (·, ·) denotes the squared geodesic disFig. 3: ML and ILS means.
tance. These parameter estimates naturally generalize their Euclidean counterparts, and they can be further shown to havemaximal likelihood when the manifold is also a symmetric space [7]. For moregeneral manifolds, like the ones under consideration in this paper, these estimates
do not attain maximal likelihood. Figure 3 show both the ILS estimate of µ andthe maximum likelihood (ML) estimate. Since the ILS estimate falls outside thesupport of the data, a significantly larger covariance matrix is needed to explain
the data, which gives a poor likelihood. To find the maximum likelihood parameters of µ and Σ we perform steepest descent directly on the data log-likelihoodusing an efficient Monte Carlo estimator of the normalization constant C [1].

3

Maximum likelihood metric learning

Determining the optimal metric (parametrized by σ) is an open question. Since
the LAND is a parametric probabilistic model it is natural to perform this model
selection using maximum likelihood. The complete data log-likelihood isN

L(σ) = −1X 2
d (xn , µ) − N log C(µ, Σ).
2 n=1 Σ(9)

It is tempting to evaluate L(σ) for several values of σ and pick the one withmaximal likelihood. This, however, is both theoretically and practically flawed.
The first issue is that the measure dM(·) used to define the LAND depends
on σ. This imply that L(σ) cannot be compared for different values of σ as they
do not rely on the same measure. The second issue is that Logµ (xn ) must be
evaluated numerically, which can become unstable when M has high curvature.
This imply that L(σ) can often not be fully evaluated when σ is small.
3.1

Likelihood bounds to cope with numerical instabilities

When numerical instabilities prevent us from evaluating L(σ) we instead relyon an easy-to-evaluate lower bound L(σ). To derive this, let vn = Logµ (xn ).
Then kvn k is the geodesic distance between µ and xn , while vn/kvn k is the initial
direction of the connecting geodesic. It is easy to provide an upper bound on thegeodesic distance by taking the length of a non-geodesic connecting curve, herechosen as the straight line connecting µ and xn . The bound then becomesZ 1p
kvn k ≤ d˜n =h(xn − µ), M(txn + (1 − t)µ)(xn − µ)idt.
(10)0

The initial orientation vn/kvn k influence the log-likelihood as the covariance Σ
is generally anisotropic. This is, however, easily bounded by picking the initial
direction as the eigenvector of Σ corresponding to the smallest eigenvalue λmin .
This then gives the final lower boundL(σ) = −N
1 X d˜2n− N log C(µ, Σ).
2 n=1 λmin

(11)

In practice, we only use the bound for data points x where the logarithm map
cannot be evaluated, and otherwise use the correct log-likelihood.
3.2

Comparing likelihoods

Since the measure dM(·) changes with σ we cannot directly compare L(σ) acrossinputs. In order to make this comparison feasible, we propose to re-normalize
the LAND with respect to the usual Lebesgue measure of the data space RD .
This amount to changing the applied measure in Eq. 6. As we lack closed-form
expressions, we perform this re-normalization using importance sampling [4]



Z
Z
exp − 21 d2Σ (µ, x)
1 2
˜C(µ, Σ) =exp − dΣ (µ, x) dx =q(x)dx (12)
2
q(x)
RD
RD


S
1X
1
1
≈
ws exp − d2Σ (µ, xs ) , xs ∼ q(x), ws =
,
(13)
S s=1
2
q(xs )
where q(x) is the proposal distribution from which we draw S samples. In our
experiments we choose q(x) = N (x|µ, Σ) with the linear mean and covarianceof the data. Thus, we ensure that the support of the proposal captures the data
manifold, but any other distribution with the desired properties can be used.

4

Results

Experimental setup: We evaluate the proposed method on both synthetic
and real data. The two-dimensional synthetic data is drawn from an arc-shaped
distribution (see Fig. 4c) [1]. We further consider features extracted from EEG
measurements during human sleep [1]; the digit “1” from MNIST; and the“Frey faces”1 . Both image modalities are projected onto their first two principal
components, and are separated into 10 and 5 folds respectively. To each data
modality, we fit a mixture of LANDs with K components.

1
Importance Sampling
Numerical Integration

0.9
0.8
0.7
0.6
0.5

unst
ab
num
ericles

0.4
0.3
0

2000

4000

6000

8000

10000

# importance samples(a) Data samples.

(c) The impact of σ.

(b) Importance sampling.
-5-10 -1

-4.7

-10 0

-10 1

-4.75

log-likelihood

log-likelihood

log-likelihood

-5.2

-5.4

-5.6

-5.8

-6

-10 2

0.07

0.09

0.11

0.13

0.15

0.17

0.5

0.19

-4.9

-4.95

-6.2
0.05

-4.8

-4.85

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

(d) Synthetic data (K=1). (e) MNIST digit 1 (K=2).

-5
0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

(f) EEG signal (K=3).

-5.2

log-likelihood

-5.4

-5.6

-5.8

-6

-6.2

-6.4
0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

(g) Frey faces (K=2).

1.5(h) Digit 1 density.

(i) EEG sleep density.

Fig. 4: Experimental results on various data sets; see text for details.

Verification: First, we validate the importance sampling scheme in Fig. 4bwhere we compare with an expensive numerical integration scheme on a predefined
grid. It is evident that importance sampling quickly gives a good approximation
to to the true normalization constant. However, choosing the correct proposal
1

http://www.cs.nyu.edu/~roweis/data.html

distribution is usually crucial for the success of the approximation [4]. Then,in Fig. 4c we show the impact of σ on the geodesic solution. When σ is small(0.01) the true geodesic cannot be computed numerically and a straight line isused to bound the likelihood (Sec. 3). For larger values of σ the geodesic canbe computed. Note that the geodesic becomes increasingly “straight” for large
values of σ.
Model selection: Figures 4d-4g show the log-likelihood bound proposed
in Sec. 3 for all data sets. In particular, we can distinguish three different regions
for the σ parameter. (1) For small values of σ the manifold has high curvatureand some geodesics cannot be computed, such that the bound penalizes the datalog-likelihood. (2) There is a range of σ values where the construction of the
manifold captures the actual underlying data structure, and in those cases weachieve the best log-likelihood. (3) For larger values of σ the manifold becomes
flat, and even if we are able compute all the geodesics the likelihood is reduced.
The reason is that when the manifolds becomes flat, significant probability mass
is assigned to regions outside of the data support, while in the other case all the
probability mass is concentrated near the data resulting to higher likelihood.

5

Discussion

Probability density estimation in non-linear spaces is essential in data analysis [6].
With the current work, we have proposed practical tools for model selection of
the metric underlying the locally adaptive normal distribution (LAND) [1]. The
basic idea amounts to picking the metric that maximize the data likelihood. A
theoretical concern is that different metrics gives different measures implying thatlikelihoods are not comparable. We have proposed to solve this by re-normalizing
according to the Lebesgue measure associated with the data space. Practically ouridea face numerical challenges when the metric has high curvature as geodesics
then become unstable to compute. Here we have proposed an easy-to-computebound on the data likelihood, which has the added benefit that metrics giving
rise to numerical instabilities are penalized. Experimental results on diverse data
sets indicate that the approach is suitable for model selection.
In this paper we have considered maximum likelihood estimation on the
training data, which can potentially overfit [10]. While we did not observe such
behavior in our experiments it is still worth investigating model selection on aheld-out test set or to put a prior on σ and pick the value that maximize theposterior probability. Both choices are straight-forward.
An interesting alternative to bounding the likelihood appears when considering
probabilistic solvers [15] for the geodesic equations (2). These represent thenumerical estimate of geodesics with a Gaussian process whose uncertainty
captures numerical approximation errors. Efficient algorithms then exist for
estimating the distribution of the geodesic arc length [3]. With these solvers,
hard-to-estimate geodesics will be associated with high variance, such that the
now-stochastic data log-likelihood also has high variance. Model selection should
then take this variance into account.

Acknowledgements. SH was supported by a research grant (15334) fromVILLUM FONDEN. LKH was supported by the Novo Nordisk Foundation
Interdisciplinary Synergy Program 2014, ’Biophysically adjusted state-informed
cortex stimulation (BASICS)’.

References
1. Arvanitidis, G., Hansen, L.K., Hauberg, S.: A locally adaptive normal distribution.
In: Advances in Neural Information Processing Systems (NIPS) (2016)
2. Belkin, M., Niyogi, P.: Laplacian eigenmaps for dimensionality reduction and data
representation. Neural computation 15(6), 1373–1396 (2003)
3. Bewsher, J., Tosi, A., Osborne, M., Roberts, S.: Distribution of Gaussian Process
Arc Lengths. In: AISTATS (2017)
4. Bishop, C.M.: Pattern Recognition and Machine Learning (Information Scienceand Statistics). Springer-Verlag New York, Inc., Secaucus, NJ, USA (2006)
5. do Carmo, M.: Riemannian Geometry. Birkhäuser (1992)
6. Chevallier, E., Barbaresco, F., Angulo, J.: Probability Density Estimation on theHyperbolic Space Applied to Radar Processing. In: GSI. pp. 753–761 (2015)
7. Fletcher, P.T.: Geodesic regression and the theory of least squares on Riemannian
manifolds. IJCV 105(2), 171–185 (2013)
8. Fletcher, P.T., Lu, C., Pizer, S.M., Joshi, S.: Principal geodesic analysis for the
study of nonlinear statistics of shape. IEEE TMI 23(8), 995–1005 (2004)
9. Freifeld, O., Hauberg, S., Black, M.J.: Model transport: Towards scalable transfer
learning on manifolds. In: CVPR (2014)
10. Hansen, L.K., Larsen, J.: Unsupervised learning and generalization. In: Neural
Networks, 1996., IEEE International Conference on. vol. 1, pp. 25–30. IEEE (1996)
11. Hauberg, S.: Principal Curves on Riemannian Manifolds. IEEE Transactions on
Pattern Analysis and Machine Intelligence (TPAMI) (2016)12. Hauberg, S., Freifeld, O., Black, M.J.: A Geometric Take on Metric Learning. In:
Advances in Neural Information Processing Systems (NIPS). pp. 2033–2041 (2012)
13. Hauberg, S., Lauze, F., Pedersen, K.S.: Unscented Kalman Filtering on Riemannian
manifolds. Journal of Mathematical Imaging and Vision 46(1), 103–120 (May 2013)
14. Hauberg, S., Pedersen, K.S.: Stick it! articulated tracking using spatial rigid object
priors. In: ACCV. LNCS, vol. 6494, pp. 758–769. Springer (2010)
15. Hennig, P., Hauberg, S.: Probabilistic Solutions to Differential Equations and their
Application to Riemannian Statistics. In: AISTATS. vol. 33 (2014)
16. Lawrence, N.: Probabilistic non-linear principal component analysis with Gaussian
process latent variable models. J. of Machine Learning Research 6, 1783–1816 (2005)
17. Pennec, X.: Intrinsic Statistics on Riemannian Manifolds: Basic Tools for GeometricMeasurements. Journal of Mathematical Imaging and Vision 25(1), 127–154 (2006)
18. Saul, L.K., Roweis, S.T.: Think globally, fit locally: unsupervised learning of low
dimensional manifolds. Journal of Machine Learning Research 4, 119–155 (2003)
19. Straub, J., Chang, J., Freifeld, O., Fisher III, J.W.: A Dirichlet Process Mixture
Model for Spherical Data. In: AISTATS (2015)
20. Tenenbaum, J.B., de Silva, V., Langford, J.C.: A Global Geometric Framework for
Nonlinear Dimensionality Reduction. Science 290(5500), 2319 (2000)
21. Tosi, A., Hauberg, S., Vellido, A., Lawrence, N.D.: Metrics for Probabilistic Geometries. In: The Conference on Uncertainty in Artificial Intelligence (UAI) (2014)
22. Zhang, M., Fletcher, P.T.: Probabilistic Principal Geodesic Analysis. In: Advancesin Neural Information Processing Systems (NIPS) 26. pp. 1178–1186 (2013)