A Locally Adaptive Normal Distribution

Georgios Arvanitidis, Lars Kai Hansen and Søren Hauberg
Technical University of Denmark, Lyngby, Denmark
DTU Compute, Section for Cognitive Systems{gear,lkai,sohau}@dtu.dk

Abstract
The multivariate normal density is a monotonic function of the distance to the mean,
and its ellipsoidal shape is due to the underlying Euclidean metric. We suggest to
replace this metric with a locally adaptive, smoothly changing (Riemannian) metricthat favors regions of high local density. The resulting locally adaptive normal
distribution (LAND) is a generalization of the normal distribution to the “manifold”
setting, where data is assumed to lie near a potentially low-dimensional manifold
embedded in RD . The LAND is parametric, depending only on a mean and acovariance, and is the maximum entropy distribution under the given metric. The
underlying metric is, however, non-parametric. We develop a maximum likelihood
algorithm to infer the distribution parameters that relies on a combination ofgradient descent and Monte Carlo integration. We further extend the LAND to
mixture models, and provide the corresponding EM algorithm. We demonstrate
the efficiency of the LAND to fit non-trivial probability distributions over bothsynthetic data, and EEG measurements of human sleep.

1

Introduction

The multivariate normal distribution is a fundamental building block in many machine learningalgorithms, and its well-known density can compactly be written as


1
2
(1)
p(x | µ, Σ) ∝ exp − distΣ (µ, x) ,2
where dist2Σ (µ, x) denotes the Mahalanobis distance for covariance matrix Σ. This distance measurecorresponds to the length of the straight line connecting µ and x, and consequently the normal
distribution is often used to model linear phenomena. When data lies near a nonlinear manifold
embedded in RD the normal distribution becomes inadequate due to its linear metric. We investigate
if a useful distribution can be constructed by replacing the linear distance function with a nonlinear
counterpart. This is similar in spirit to Isomap [21] that famously replace the linear distance with ageodesic distance measured over a neighborhood graph spanned by the data, thereby allowing fora nonlinear model. This is, however, a discrete distance measure that is only well-defined over the
training data. For a generative model, we need a continuously defined metric over the entire RD .
Following Hauberg et al. [9] we learn a smoothly changing metric that favors regions of high density
i.e., geodesics tend to move near the data. Under this metric, the data space is interpreted as a
D-dimensional Riemannian manifold. This “manifold learning” does not change dimensionality, but
merely provides a local description of the data. The Riemannian view-point, however, gives a strong
mathematical foundation upon which the proposed distribution can be developed. Our work, thus,
bridges work on statistics on Riemannian manifolds [15, 23] with manifold learning [21].
We develop a locally adaptive normal distribution (LAND) as follows: First, we construct a metricthat captures the nonlinear structure of the data and enables us to compute geodesics; from this, an
30th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain.

Geodesics

LAND model

Data
LAND mean
LinearGeodesic

Linear model

LAND mean

Linear mean

Figure 1: Illustration of the LAND using MNIST images of the digit 1 projected onto the first 2
principal components. Left: comparison of the geodesic and the linear distance. Center: the proposed
locally adaptive normal distribution. Right: the Euclidean normal distribution.
unnormalized density is trivially defined. Second, we propose a scalable Monte Carlo integration
scheme for normalizing the density with respect to the measure induced by the metric. Third, wedevelop a gradient-based algorithm for maximum likelihood estimation on the learned manifold. We
further consider a mixture of LANDs and provide the corresponding EM algorithm. The usefulnessof the model is verified on both synthetic data and EEG measurements of human sleep stages.
Notation: all points x ∈ RD are considered as column vectors, and they are denoted with boldDlowercase characters. S++
represents the set of symmetric D × D positive definite matrices. The
learned Riemannian manifold is denoted M, and its tangent space at x ∈ M is denoted Tx M.

2

A Brief Summary of Riemannian Geometry

We start our exposition with a brief review of Riemannian manifolds [6]. These smooth manifolds arenaturally equipped with a distance measure, and are commonly used to model physical phenomena
such as dynamical or periodic systems, and many problems that have a smooth behavior.
D
Definition 1. A smooth manifold M together with a Riemannian metric M : M → S++
is calleda Riemannian manifold. The Riemannian metric M encodes a smoothly changing inner producthu, M(x)vi on the tangent space u, v ∈ Tx M of each point x ∈ M.

Remark 1. The Riemannian metric M(x) acts on tangent vectors, and may, thus, be interpreted asa standard Mahalanobis metric restricted to an infinitesimal region around x.
The local inner product based on M is a suitable model for capturing local behavior of data, i.e.
manifold learning. From the inner product, we can define geodesics as length-minimizing curvesconnecting two points x, y ∈ M, i.e.
Z 1pγ̂ = argmin
hγ 0 (t), M(γ(t))γ 0 (t)idt, s.t. γ(0) = x, γ(1) = y.
(2)
γ

0

Here M(γ(t)) is the metric tensor at γ(t), and the tangent vector γ 0 denotes the derivative (velocity) of γ. The distance between x and y is defined as the length of thegeodesic. A standard result from differential geometry is that the geodesic can be foundas the solution to a system of 2nd order ordinary differential equations (ODEs) [6, 9]:|1∂vec[M(γ(t))]γ (t) = − M−1 (γ(t))(γ 0 (t) ⊗ γ 0 (t))
2
∂γ(t)
(3)

v = Logx (y)

γ(t)

00

x

subject to γ(0) = x, γ(1) = y. Here vec[·] stacks the columnsof a matrix into a vector and ⊗ is the Kronecker product.
y = Expx (v)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 Figure 2: An illustration of the exγ(t) = Expx (t · v) is a geodesic originating at x with initial ponential and logarithmic maps.
2

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 fromx to y. These operations are illustrated in Fig. 2. The exponential and the logarithmic map can
be computed by solving Eq. 3 numerically, as an initial value problem (IVP) or a boundary value
problem (BVP) respectively. In practice the IVPs are substantially faster to compute than the BVPs.
The Mahalanobis distance is naturally extended to Riemannian manifolds as dist2Σ (x, y) =hLogx (y), Σ−1 Logx (y)i. From this, Pennec [15] considered the Riemannian normal distribution


1
1
−1
pM (x | µ, Σ) = exp − hLogµ (x), Σ Logµ (x)i , x ∈ M(4)C2
and showed that it is the manifold-valued distribution with maximum entropy subject to a known
mean and covariance. This distribution is an instance of Eq. 1 and is the distribution we consider inthis paper. Next, we consider standard “intrinsic least squares” estimates of µ and Σ.
2.1

Intrinsic Least Squares Estimators

Let the data be generated from an unknown probability distribution qM (x) on a manifold. Then it iscommon [15] to define the intrinsic mean of the distribution as the point that minimize the varianceZµ̂ = argmin
dist2 (µ, x)qM (x)dM(x),
(5)
µ∈M

M

where dM(x) is the measure (or infinitesimal volume element) induced by the metric. Based on themean, a covariance matrix can be definedZΣ̂ =Logµ̂ (x)Logµ̂ (x)| qM (x)dM(x),
(6)
D(µ̂)

where D(µ̂) is the domain over which Tµ̂ M is well-defined. For the manifolds we consider, thedomain D(µ̂) is RD . Practical estimators of µ̂ rely on gradient-based optimization to find a local
minimizer of Eq. 5, which is well-defined [12]. For finite data {xn }N
n=1 , the descent direction is
PN
proportional to v̂ = n=1 Logµ (xn ) ∈ Tµ M, and the updated mean is a point on the geodesic
curve γ(t) = Expµ (t · v̂). After estimating the mean, the empirical covariance matrix is estimatedPN
as Σ̂ = N 1−1 n=1 Logµ̂ (xn )Logµ̂ (xn )| . It is worth noting that even though these estimators arenatural, they are not maximum likelihood estimates for the Riemannian normal distribution (4).
In practice, the intrinsic mean often falls in regions of low data density [8]. For instance, consider
data distributed uniformly on the equator of a sphere, then the optima of Eq. 5 is either of the poles.
Consequently, the empirical covariance is often overestimated.

3

A Locally Adaptive Normal Distribution

We now have the tools to define a locally adaptive normal distribution (LAND): we replace thelinear Euclidean distance with a locally adaptive Riemannian distance and study the corresponding
Riemannian normal distribution (4). By learning a Riemannian manifold and using its structure toestimate distributions of the data, we provide a new and useful link between Riemannian statistics
and manifold learning.
3.1

Constructing a MetricIn the context of manifold learning, Hauberg et al. [9] suggest to model the local behavior of the data
manifold via a locally-defined Riemannian metric. Here we propose to use a local covariance matrix
to represent the local structure of the data. We only consider diagonal covariances for computational
efficiency and to prevent the overfitting. The locality of the covariance is defined via an isotropic
Gaussian kernel of size σ. Thus, the metric tensor at x ∈ M is defined as the inverse of a local
diagonal covariance matrix with entries
!−1
!
N
2
X
kxn − xk2
2
Mdd (x) =wn (x)(xnd − xd ) + ρ, with wn (x) = exp −
. (7)
2σ 2
n=1
3

Here xnd is the dth dimension of the nth observation, and ρ a regularization parameter to avoid
singular covariances. This defines a smoothly changing (hence Riemannian) metric that captures the
local structure of the data. It is easy to see that if x is outside of the support of the data, then themetric tensor is large. Thus, geodesics are “pulled” towards the data where the metric is small. Notethat the proposed metric is not invariant to linear transformations.While we restrict our attention tothis particular choice, other learned metrics are equally applicable, c.f. [22, 9].
3.2

Estimating the Normalization Constant

The normalization constant of Eq. 4 is by definitionZ
1
C(µ, Σ) =exp − hLogµ (x), Σ−1 Logµ (x)i dM(x),
2
M

(8)

where dM(x) denotes the measure induced by the Riemannian metric. The constant C(µ, Σ) depends
not only on the covariance matrix, but also on the mean of the distribution, and the curvature of the
manifold (captured by the logarithm map). For a general learned manifold, C(µ, Σ) is inaccessible inclosed-form and we resort to numerical techniques. We start by rewriting Eq. 8 as


Z
q
1
M(Expµ (v)) exp − hv, Σ−1 vi dv.
C(µ, Σ) =
(9)
2
Tµ M
In effect, we integrate the distribution over the tangent space Tµ M instead of directly over the
manifold. This transformation relies on the fact that the volume of an infinitely small area on
the manifold can be computed in the tangent space if we take the deformation of the metric into
account[15]. This deformation is captured by the measure which, in the tangent space,qq is dM(x) =M(Expµ (v)) dv. For notational simplicity we define the function m(µ, v) =M(Expµ (v)) ,which intuitively captures the cost for a point to be outside the data support (m is large in low density
areas and small where the density is high).
We estimate the normalization constant (9) using Monte Carlo integration. We first multiply and divide the integral with thep normalization constant of the Euclidean normal distribution Z = (2π)D |Σ|.
Then, the integral becomes an expectation estimation problem
C(µ, Σ) = Z · EN (0,Σ) [m(µ, v)], which can be estimated numerically as

Intrinsic
Least
SquaresS

C(µ, Σ) 'ZX
m(µ, vs ),S s=1

where vs ∼ N (0, Σ)

(10)

LAND

Figure 3: Comparison ofand S is the number of samples on Tµ M. The computationally LAND and intrinsic least
expensive element is to evaluate m, which in turn requires evaluating squares means.
Expµ (v). This amounts to solving an IVP numerically, which isfairly fast. Had we performed the integration directly on the manifold (8) we would have had toevaluate the logarithm map, which is a much more expensive BVP. The tangent space integration,
thus, scales better.
3.3

Inferring Parameters

Assuming an independent and identically distributed dataset {xn }N
n=1 , we can write their joint
QN
distribution as pM (x1 , . . . , xN ) = n=1 pM (xn | µ, Σ). We find parameters µ and Σ by maximum
likelihood, which we implement by minimizing the mean negative log-likelihoodN

1 X{µ̂, Σ̂} = argmin φ (µ, Σ) = argmin
hLogµ (xn ), Σ−1 Logµ (xn )i + log (C(µ, Σ)) .
µ∈Mµ∈M 2N n=1
D
Σ∈S++

D
Σ∈S++

(11)
D
The first term of the objective function φ : M × S++
is a data-fitting term, while the second can beseen as a force that both pulls the mean closer to the high density areas and shrinks the covariance.
Specifically, when the mean is in low density areas, as well as when the covariance gives significant
4probability to those areas, the value of m(µ, v)
Algorithm 1 LAND maximum likelihood
will by construction be large. Consequently,
C(µ, Σ) will increase and these solutions will be Input: the data {xn }N
n=1 , stepsizes αµ , αA
penalized. In practice, we find that the maximum Output: the estimated µ̂, Σ̂, C(ˆ µ̂, Σ̂)
likelihood LAND mean generally avoids low den- 1: initialize µ0 , Σ0 and t ← 0sity regions, which is in contrast to the standard
2: repeat
intrinsic least squares mean (5), see Fig. 3.
3:estimate C(µt , Σt ) using Eq. 10
compute dµ φ(µt , Σt ) using Eq. 12
In practice we optimize φ using block coordinate 4:
descent: we optimize the mean keeping the co- 5:µt+1 ← Expµt (αµ dµ φ(µt , Σt ))variance fixed and vice versa. Unfortunately, both 6:
estimate C(µt+1 , Σt ) using Eq. 10of the sub-problems are non-convex, and unlike
7:
compute ∇A φ(µt+1 , Σt ) using Eq. 13
the linear normal distribution, they lack a closedAt+1 ← At − αA ∇A φ(µt+1 , Σt )form solution. Since the logarithm map is a dif- 8:
9:
Σt+1 ← [(At+1 )| At+1 ]−1ferentiable function, we can use gradient-basedt←t+1
techniques to infer µ and Σ. Below we give the 10:
t+1
t 2
descent direction for µ and Σ and the correspond- 11: until φ(µt+1 , Σ ) − φ(µt , Σ ) 2 ≤ 
ing optimization scheme is given in Algorithm 1.
Initialization is discussed in the supplements.
Optimizing µ: the objective function is differentiable with respect to µ [6], and using that
−1∂
Logµ (x)i = −2Σ−1 Logµ (x), we get the gradient∂µ hLogµ (x), Σ
#
"
N
S
X
Z
1 X
−1
Logµ (xn ) −m(µ, vs )vs .
(12)∇µ φ(µ, Σ) = −ΣN n=1
C(µ, Σ) · S s=1
It is easy to see that this gradient is highly dependent on the condition number of Σ. We find that this,at times, makes the gradient unstable, and choose to use the steepest descent direction instead of thegradient direction. This is equal to dµ φ(µ, Σ) = −Σ∇µ φ(µ, Σ) (see supplements).
D
Optimizing Σ: since the covariance matrix by definition is constrained to be in the space S++,a−1|common trick is to decompose the matrix as Σ = A A, and optimize the objective with respect
to A. The gradient of this factor is (see supplements for derivation)
"
#
N
S
X
Z
1 X
|
|
Logµ (xn )Logµ (xn ) −m(µ, vs )vs vs .
(13)
∇A φ(µ, Σ) = A
N n=1
C(µ, Σ) · S s=1

Here the first term fits the given data by increasing the size of the covariance matrix, while the second
term regularizes the covariance towards a small matrix.
3.4

Mixture of LANDs

At this point we can find maximum likelihood estimates of the LAND model. We can easily extend
this to mixtures of LANDs: Following the derivation of the standard Gaussian mixture model [3], ourobjective function for inferring the parameters of the LAND mixture model is formulated as follows


K X
N
X
1
−1
hLogµk (xn ), Σk Logµk (xn )i + log(C(µk , Σk )) − log(πk ) , (14)
ψ(Θ) =
rnk
2
n=1
k=1

π p (x | µk ,Σk )PKk M nwhere Θ = {µk , Σk }Kis the probability that xn is generated by the
k=1 , rnk =l=1 πl pM (xn | µl ,Σl )PKk th component, and k=1 πk = 1, πk ≥ 0. The corresponding EM algorithm is in the supplements.

4

Experiments

In this section we present both synthetic and real experiments to demonstrate the advantages of the
LAND. We compare our model with both the Gaussian mixture model (GMM), and a mixture of
LANDs using least squares (LS) estimators (5, 6). Since the latter are not maximum likelihood
estimates we use a Riemannian K-means algorithm to find cluster centers. In all experiments we
use S = 3000 samples in the Monte Carlo integration. This choice is investigated empirically in the
supplements. Furthermore, we choose σ as small as possible, while ensuring that the manifold issmooth enough that geodesics can be computed numerically.
54.1

Synthetic Data Experiments

Mean negative log-likelihood

7
As a first experiment, we generate a nonlinear data-manifold
GMMLS
by sampling from a mixture of 20 Gaussians positioned along a
6
LAND
True
half-ellipsoidal curve (see left panel of Fig. 5). We generate 10
5datasets with 300 points each, and fit for each dataset the three4
models with K = 1, . . . , 4 number of components. Then, we
3
generate 10000 samples from each fitted model, and we compute the mean negative log-likelihood of the true generative2
distribution using these samples. Fig. 4 shows that the LAND
1
learns faster the underlying true distribution, than the GMM.
0
Moreover, the LAND perform better than the least squares esti12
3
4
Number of mixture components
mators, which overestimates the covariance. In the supplements
we show, using the standard AIC and BIC criteria, that the op- Figure 4: The mean negative logtimal LAND is achieved for K = 1, while for the least squares likelihood experiment.
estimators and the GMM, the optimal is achieved for K = 3
and K = 4 respectively.

In addition, in Fig. 5 we show the contours for the LAND and the GMM for K = 2. There, we
can observe that indeed, the LAND adapts locally to the data and reveals their underlying nonlinear
structure. This is particularly evident near the “boundaries” of the data-manifold.
Geodesics

Gaussian mixture model

LAND mixture model
LAND mean

GMM mean

Data
LAND means
Geodesics, cluster 1
Geodesics, cluster 2

Figure 5: Synthetic data and the fitted models. Left: the given data, the intensity of the geodesics
represent the responsibility of the point to the corresponding cluster. Center: the contours of the
LAND mixture model. Right: the contours of the Gaussian mixture model.
We extend this experiment to a clustering task (see left panel of Fig. 6 for data). The center and rightpanels of Fig. 6 show the contours of the LAND and Gaussian mixtures, and it is evident that the
LAND is substantially better at capturing non-ellipsoidal clusters. Due to space limitations, we movefurther illustrative experiments to the supplementary material and continue with real data.
4.2

Modeling Sleep Stages

We consider electro-encephalography (EEG) measurements of human sleep from 10 subjects, part of
the PhysioNet database [11, 7, 5]. For each subject we get EEG measurements during sleep fromtwo electrodes on the front and the back of the head, respectively. Measurements are sampled at
fs = 100Hz, and for each 30 second window a so-called sleep stage label is assigned from the set{1, 2, 3, 4, REM, awake}. Rapid eye movement (REM) sleep is particularly interesting, characterizedby having EEG patterns similar to the awake state but with a complex physiological pattern, involvinge.g., reduced muscle tone, rolling eye movements and erection [16]. Recent evidence points to theimportance of REM sleep for memory consolidation [4]. Periods in which the sleeper is awake are
typically happening in or near REM intervals. Thus we here consider the characterization of sleep interms of three categories REM, awake, and non-REM, the latter a merger of sleep stages 1 − 4.
We extract features from EEG measurements as follows: for each subject we subdivide the 30 second
windows to 10 seconds, and apply a short-time-Fourier-transform to the EEG signal of the frontal
electrode with 50% overlapping windows. From this we compute the log magnitude of the spectrumlog(1 + |f |) of each window. The resulting data matrix is decomposed using Non-Negative Matrix
Factorization (10 random starts) into five factors, and we use the coefficients as 5D features. In Fig. 7
we illustrate the nonlinear manifold structure based on a three factor analysis.
6Geodesics

LAND mixture model

Gaussian mixture model

LAND mean

GMM mean

Data
LAND means
Geodesics, cluster 1
Geodesics, cluster 2
Geodesics

LAND mixture model

Gaussian mixture model

LAND mean

GMM mean

Data
LAND means
Geodesics, cluster 1
Geodesics, cluster 2

Figure 6: The clustering problem for two synthetic datasets. Left: the given data, the intensity of thegeodesics represent the responsibility of the point to the corresponding cluster. Center: the LAND
mixture model. Right: the Gaussian mixture model.

We perform clustering on the data and evaluate the alignment
between cluster labels and sleep stages using the F-measure[14]. The LAND depends on the parameter σ to construct the
1-4
R.E.M.
metric tensor, and in this experiment it is less straightforward to
awake
select σ because of significant intersubject variability. First, we
fixed σ = 1 for all the subjects. From the results in Table 1 we
observe that for σ = 1 the LAND(1) generally outperforms the
GMM and achieves much better alignment. To further illustrate
the effect of σ we fitted a LAND for σ = [0.5, 0.6, . . . , 1.5]
and present the best result achieved by the LAND. Selecting σ
this way leads indeed to higher degrees of alignment further unFigure 7: The 3 leading factors forderlining that the conspicuous manifold structure and the rather
subject “s151”.
compact sleep stage distributions in Fig. 7 are both captured
better with the LAND representation than with a linear GMM.
Table 1: The F-measure result for 10 subjects (the closer to 1 the better).

5

s001

s011

s042

s062

s081

s141

s151

s161

s162

s191

LAND(1)
GMM

0.831
0.812

0.701
0.690

0.670
0.675

0.740
0.651

0.804
0.798

0.870
0.870

0.820
0.794

0.780
0.775

0.747
0.747

0.786
0.776

LAND

0.831

0.716

0.695

0.740

0.818

0.874

0.830

0.783

0.750

0.787

Related WorkWe are not the first to consider Riemannian normal distributions, e.g. Pennec [15] gives a theoreticalanalysis of the distribution, and Zhang and Fletcher [23] consider the Riemannian counterpart ofprobabilistic PCA. Both consider the scenario where the manifold is known a priori. We adapt the
distribution to the “manifold learning” setting by constructing a Riemannian metric that adapts to thedata. This is our overarching contribution.
Traditionally, manifold learning is seen as an embedding problem where a low-dimensional representation of the data is sought. This is useful for visualization [21, 17, 18, 1], clustering [13],semi-supervised learning [2] and more. However, in embedding approaches, the relation between a7new point and the embedded points are less well-defined, and consequently these approaches are lesssuited for building generative models. In contrast, the Riemannian approach gives the ability to measure continuous geodesics that follow the structure of the data. This makes the learned Riemannian
manifold a suitable space for a generative model.
Simo-Serra et al. [19] consider mixtures of Riemannian normal distributions on manifolds thatare known a priori. Structurally, their EM algorithm is similar to ours, but they do not accountfor the normalization constants for different mixture components. Consequently, their approach isinconsistent with the probabilistic formulation. Straub et al. [20] consider data on spherical manifolds,
and further consider a Dirichlet process prior for determining the number of components. Such aprior could also be incorporated in our model. The key difference to our work is that we consider
learned manifolds as well as the following complications.

6

Discussion

In this paper we have introduced a parametric locally adaptive normal distribution. The idea is toreplace the Euclidean distance in the ordinary normal distribution with a locally adaptive nonlinear
distance measure. In principle, we learn a non-parametric metric space, by constructing a smoothly
changing metric that induces a Riemannian manifold, where we build our model. As such, we proposea parametric model over a non-parametric space.
The non-parametric space is constructed using a local metric that is the inverse of a local covariancematrix. Here locality is defined via a Gaussian kernel, such that the manifold learning can be seen
as a form of kernel smoothing. This indicates that our scheme for learning a manifold might notscale to high-dimensional input spaces. In these cases it may be more practical to learn the manifold
probabilistically [22] or as a mixture of metrics [9]. This is feasible as the LAND estimation procedureis agnostic to the details of the learned manifold as long as exponential and logarithm maps can beevaluated.
Once a manifold is learned, the LAND is simply a Riemannian normal distribution. This is a naturalmodel, but more intriguing, it is a theoretical interesting model since it is the maximum entropy
distribution for a fixed mean and covariance [15]. It is generally difficult to build locally adaptive
distributions with maximum entropy properties, yet the LAND does this in a fairly straight-forwardmanner. This is, however, only a partial truth as the distribution depends on the non-parametric space.
The natural question, to which we currently do not have an answer, is whether a suitable maximum
entropy manifold exist?
Algorithmically, we have proposed a maximum likelihood estimation scheme for the LAND. This
combines a gradient-based optimization with a scalable Monte Carlo integration method. Once
exponential and logarithm maps are available, this procedure is surprisingly simple to implement. We
have demonstrated the algorithm on both real and synthetic data and results are encouraging. We
almost always improve upon a standard Gaussian mixture model as the LAND is better at capturing
the local properties of the data.
We note that both the manifold learning aspect and the algorithmic aspect of our work can be improved.
It would be of great value to learn the parameter σ used for smoothing the Riemannian metric, and ingeneral, more adaptive learning schemes are of interest. Computationally, the bottleneck of our work
is evaluating the logarithm maps. This may be improved by specialized solvers, e.g. probabilistic
solvers [10], or manifold-specific heuristics.
The ordinary normal distribution is a key element in many machine learning algorithms. We expectthat many fundamental generative models can be extended to the “manifold” setting simply by
replacing the normal distribution with a LAND. Examples of this idea include Naïve Bayes, Linear
Discriminant Analysis, Principal Component Analysis and more. Finally we note that standard
hypothesis tests also extend to Riemannian normal distributions [15] and hence also to the LAND.
Acknowledgements. LKH was funded in part by the Novo Nordisk Foundation Interdisciplinary
Synergy Program 2014, ’Biophysically adjusted state-informed cortex stimulation (BASICS)’. SH
was funded in part by the Danish Council for Independent Research, Natural Sciences.

8

References[1] M. Belkin and P. Niyogi. Laplacian Eigenmaps for Dimensionality Reduction and Data Representation.
Neural Computation, 15(6):1373–1396, June 2003.
[2] M. Belkin, P. Niyogi, and V. Sindhwani. Manifold Regularization: A Geometric Framework for Learningfrom Labeled and Unlabeled Examples. JMLR, 7:2399–2434, Dec. 2006.
[3] C. M. Bishop. Pattern Recognition and Machine Learning (Information Science and Statistics). SpringerVerlag New York, Inc., Secaucus, NJ, USA, 2006.
[4] R. Boyce, S. D. Glasgow, S. Williams, and A. Adamantidis. Causal evidence for the role of REM sleep
theta rhythm in contextual memory consolidation. Science, 352(6287):812–816, 2016.
[5] A. Delorme and S. Makeig. EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics
including independent component analysis. J. Neurosci. Methods, page 21, 2004.
[6] M. do Carmo. Riemannian Geometry. Mathematics (Boston, Mass.). Birkhäuser, 1992.
[7] A. L. Goldberger, L. A. N. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B.
Moody, C.-K. Peng, and H. E. Stanley. PhysioBank, PhysioToolkit, and PhysioNet: Components of a New
Research Resource for Complex Physiologic Signals. Circulation, 101(23):e215–e220, 2000 (June 13).
[8] S. Hauberg. Principal Curves on Riemannian Manifolds. IEEE Transactions on Pattern Analysis andMachine Intelligence (TPAMI), 2016.
[9] S. Hauberg, O. Freifeld, and M. J. Black. A Geometric Take on Metric Learning. In Advances in Neural
Information Processing Systems (NIPS) 25, pages 2033–2041, 2012.
[10] P. Hennig and S. Hauberg. Probabilistic Solutions to Differential Equations and their Application to
Riemannian Statistics. In Proceedings of the 17th international Conference on Artificial Intelligence andStatistics (AISTATS), volume 33, 2014.
[11] S. A. Imtiaz and E. Rodriguez-Villegas. An open-source toolbox for standardized use of PhysioNet Sleep
EDF Expanded Database. In 2015 37th Annual International Conference of the IEEE Engineering in
Medicine and Biology Society (EMBC), pages 6014–6017, Aug 2015.
[12] H. Karcher. Riemannian center of mass and mollifier smoothing. Communications on Pure and Applied
Mathematics, 30(5):509–541, 1977.
[13] U. Luxburg. A Tutorial on Spectral Clustering. Statistics and Computing, 17(4):395–416, Dec. 2007.
[14] R. Marxer, H. Purwins, and A. Hazan. An f-measure for evaluation of unsupervised clustering with
non-determined number of clusters. Report of the EmCAP project (European Commission FP6-IST), pages
1–3, 2008.
[15] X. Pennec. Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements.
Journal of Mathematical Imaging and Vision, 25(1):127–154, July 2006.
[16] D. Purves, G. Augustine, D. Fitzpatrick, W. Hall, A. LaMantia, J. McNamara, and L. White. Neuroscience,
2008. De Boeck, Sinauer, Sunderland, Mass.
[17] S. T. Roweis and L. K. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science,
290:2323–2326, 2000.
[18] 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.
[19] E. Simo-Serra, C. Torras, and F. Moreno-Noguer. Geodesic Finite Mixture Models. In Proceedings of the
British Machine Vision Conference. BMVA Press, 2014.
[20] J. Straub, J. Chang, O. Freifeld, and J. W. Fisher III. A Dirichlet Process Mixture Model for Spherical
Data. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2015.
[21] J. B. Tenenbaum, V. de Silva, and J. C. Langford. A Global Geometric Framework for Nonlinear
Dimensionality Reduction. Science, 290(5500):2319, 2000.
[22] 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.
[23] M. Zhang and P. Fletcher. Probabilistic Principal Geodesic Analysis. In Advances in Neural Information
Processing Systems (NIPS) 26, pages 1178–1186, 2013.

9