Harmonic Exponential Families on Manifolds

Taco S. Cohen
University of Amsterdam

T. S . COHEN @ UVA . NL

Max Welling
University of Amsterdam
University of California Irvine
Canadian Institute for Advanced Research

M . WELLING @ UVA . NL

Abstract
In a range of fields including the geosciences,
molecular biology, robotics and computer vision,
one encounters problems that involve randomvariables on manifolds. Currently, there is a lackof flexible probabilistic models on manifolds thatare fast and easy to train. We define an extremely
flexible class of exponential family distributions
on manifolds such as the torus, sphere, and rotation groups, and show that for these distributions the gradient of the log-likelihood can becomputed efficiently using a non-commutative
generalization of the Fast Fourier Transform(FFT). We discuss applications to Bayesian camera motion estimation (where harmonic exponential families serve as conjugate priors), and modelling of the spatial distribution of earthquakes
on the surface of the earth. Our experimental results show that harmonic densities yield a significantly higher likelihood than the best competing
method, while being orders of magnitude fasterto train.

1. Introduction
Many problems in science and engineering involve random variables on manifolds. In the geosciences, for example, one deals with measurements such as the locationsof earthquakes and weather events on the spherical surfaceof the earth. In robotics and computer vision, unobserved
Lie group transformations such as rotations and rigid-body
motions play an important role in motion understanding,
localization and alignment problems. Classical probabilistic models cannot be applied to data on manifolds because
Proceedings of the 32 nd International Conference on Machine
Learning, Lille, France, 2015. JMLR: W&CP volume 37. Copyright 2015 by the author(s).

these models do not respect the manifold topology (having discontinuities at the value 0 ≡ 2π of a circular variable, for instance), and because they are not equivariant (amanifold-preserving transformation of the data would takethe distribution outside the model family). Among manifold distributions there are currently none that are both flexible and efficiently trainable.
In this paper we study a very flexible class of densities
on compact Lie groups (such as the group of rotations in
two or three dimensions) and homogeneous spaces (suchas the circle, torus, and sphere). We refer to these densities as harmonic exponential families, because they arebased on a generalized form of Fourier analysis known
as non-commutative harmonic analysis (Chirikjian & Kyatkin, 2001). Specifically, the sufficient statistics of thesefamilies are given by functions that are analogous to thesinusoids on the circle.
We show that the moment map (Wainwright & Jordan,
2007) that takes the natural parameters of an exponential
family and produces the moments of the distribution can becomputed very efficiently for harmonic exponential families using generalized Fast Fourier Transform (FFT) algorithms. This leads directly to a very efficient maximumlikelihood estimation procedure that is applicable to any
manifold for which an FFT algorithm has been developed,
and that enjoys the convexity and convergence propertiesof exponential families.
We apply the harmonic exponential family of the sphere to
the problem of modeling the spatial distribution of earthquakes on the surface of the earth. Our model significantly
outperforms the best competing method (a mixture model),while being orders of magnitude faster to train.
Harmonic exponential families also arise naturally as conjugate priors in a Bayesian framework for estimating transformations from correspondence pairs. In this case thepoints on the manifold (the transformations) are not observed directly. Instead we see a pair of vectors x and yHarmonic Exponential Families— images, point clouds, or other data — that provide information about the transformation that produced one fromthe other. If we have a prior p(g) over transformations anda likelihood p(y | x, g) that measures how likely y is to be
the g-transformed version of x, we can consider the posterior over transformations p(g | x, y). Typically, the posterior turns into a complicated and intractable distribution,
but we show that if the prior is a harmonic density and thelikelihood is Gaussian, the posterior distribution is again aharmonic density whose parameters are easily obtained using the generalized FFT algorithm. Furthermore, the global
mode of this posterior (the optimal transformation) can becomputed efficiently by performing yet another FFT.
In this paper we provide both an abstract treatment of thetheory of harmonic densities that covers a fairly broad classof manifolds in a uniform and easily understandable manner, and a concrete instantiation of this theory in the caseof the rotation groups SO(2) and SO(3), and their homogeneous spaces, the circle S 1 and sphere S 2 .
The rest of this paper is organized as follows. We begin by discussing related work in section 2, followed bya brief summary of non-commutative harmonic analysis
for a machine learning audience in section 3. In section
4, we define harmonic exponential families and present an
FFT-based maximum-likelihood estimation algorithm. In
section 5, we show that harmonic exponential families are
the conjugate priors in the Bayesian transformation estimation problem, and present an efficient MAP inference
algorithm. Our earthquake modelling experiments are presented in section 6, followed by a discussion and conclusion.

2. Related Work
The harmonic exponential families were first defined abstractly (but not named) by Diaconis (1988). Despite theirlong history and a significant body of literature devoted tothem (see Mardia & Jupp (1999)), this class of densities
has remained intractable until now for all but the simplest
cases.
In fact, various commonly used distributions on the circle and the sphere are harmonic densities of low degree.
Among these are the 2-parameter von-Mises distribution
on the circle and the 5-parameter Kent distribution on thesphere. These are both exponential families, about which
Mardia & Jupp (1999) write (p. 175): “Although exponential models have many pleasant inferential properties,
the need to evaluate the normalizing constant (or at least
the first derivative of its logarithm) can be a practical difficulty.”
What is a practical difficulty for the unimodal distributions
with few parameters mentioned above becomes a show-

stopper for more flexible exponential families with many
parameters. The harmonic exponential family for the circle
is known as the generalized von-Mises distribution, and canbe defined for any band-limit / order / degree L. However,no scalable maximum-likelihood estimation algorithm isknown. Gatto & Jammalamadaka (2007) (who work onlywith the 4-parameter L = 2 distribution) compute the normalizing constant by a truncated infinite sum, where eachterm involves an expensive Bessel function evaluation.
Beran (1979) studies an exponential family that is equivalent to the harmonic exponential family on the sphere, butdoes not provide a scalable learning algorithm. At each
iteration, the proposed algorithm computes all O(L2 ) moments of the distribution by numerical integration. Thisrequires O(L2 ) samples per integration, making the periteration cost O(L4 ). Beran further suggests using a second
order optimization method, which would further increase
the per-iteration cost to O(L6 ).
This is clearly not feasible when L is measured in the hundreds, and parameter counts in the 10’s of thousands, as isneeded in the experiments reported in section 6. The algorithm described in this paper is simple, generic across
manifolds, and fast (per-iteration complexity O(L2 log2 L)
in the spherical case, for instance). It can be applied to any
manifold for which an FFT has been developed.

3. Preliminaries
The manifolds we consider in this paper are either Liegroups or closely related manifolds called homogeneous
spaces, and the sufficient statistics that we use come fromharmonic analysis on these manifolds. Since these concepts are not widely known in machine learning, we will
review them in this section. For more details, we refer
the reader to Chirikjian & Kyatkin (2001); Sugiura (1990);
Kondor (2008); Goodman & Wallach (2009).
3.1. Lie Groups
A transformation group G is a set of invertible transformations that is closed under composition and taking inverses:
for any g, h ∈ G, the composition gh is again a member ofG, and so are the inverses g −1 and h−1 .
A Lie group is a group that is also a differentiable manifold. For example, the group SO(3) of 3D rotations is a
Lie group. It can be represented as a set of matrices,SO(3) = {R ∈ R3×3 | RRT = I, det(R) = 1},(1)

so one way to think about SO(3) is as a 3-dimensional manifold embedded in R3×3 . For technical reasons, we willfurther restrict our attention to compact Lie groups, i.e. Liegroups that are closed and bounded.

Harmonic Exponential Families

3.2. Harmonic analysis on compact Lie groups
The basic idea of the generalized Fourier transform oncompact Lie groups is to expand a function f : G → C
as a linear combination of carefully chosen basis functions.
These basis functions have very special properties, because
they are the matrix elements of irreducible unitary representations (IURs) of G. These terms will now be explained.
A representation of a group G on a vector space V is a
map R from the group to the set of invertible linear transformations of V that preserves the group structure in thefollowing sense:R(gh) = R(g)R(h).

(2)

Note that gh denotes composition of group elements while
R(g)R(h) denotes matrix multiplication (once we choosea basis for V , that is).

that standard Fourier analysis is a special case of the moregeneral transform to be defined shortly.
The matrix elements of IURs of SO(3) are known as
Wigner D-functions. They are defined for λ = 0, 1, 2, . . .
and −λ ≤ m, n ≤ λ, so the irreducible representationsare (2λ + 1)-dimensional. Wigner D-functions can be expressed as sums over products of sinusoids or complex exponentials (Pinchon & Hoggan, 2007), but the formulae aresomewhat unwieldy so that it is easier to think only about
their general properties.
The most important general property of the matrix elementsof IURs is that they are orthogonal:Zλλ0λλ0 (g) dµ(g)hUmn(g), Um
Umn(g) Um0 n0 (g)i ≡
0 n0
G
(4)
δλλ0 δmm0 δnn0
=
.
dim λ

In computer vision, we encounter group representations in
the following way. An image is represented as a vectorx ∈ Rn of pixel intensities, and to be concrete, we take G
to be the group SO(2) consisting of rotations of the plane.
Then R(θ) is the matrix such that R(θ)x is the image x
rotated by angle θ.

Here µ is the normalized Haar measure, which is the natural
way to measure volumes in G (Sugiura, 1990), and dim λis the dimension of the representation. One can verify thatthe complex exponentials eiλg are indeed orthonormal withdg
and dim λ = 1.
dµ(g) = 2π

As this example shows, many representations do notchange the norm of the vectors on which they act: ∀g ∈G, ∀x ∈ V : kU (g) xk = kxk. Such representations arecalled unitary. Unitary representations tend to be easier towork with both analytically and numerically.

Intuitively, the matrix elements are like a “complete orthogonal basis” for the space L2 (G) of square integrable
functions on G. That is, it can be proven that any functionf ∈ L2 (G) can be written as
XX


λ
λ
f (g) =ηmn
Tmn(g) ≡ F −1 η (g),(5)

Given a unitary representation U and a unitary matrix
F , one can define an equivalent representation T (g) =F −1 U (g)F . In computer vision one can think of F as amatrix containing image features in the rows. One can now
try to find F such that for every g, the matrix F −1 U (g)Fis block diagonal with the same block structure. If we continue to block-diagonalize until no further diagonalizationis possible, we end up with blocks called irreducible unitary representations1 . The IURs of a compact group can beindexed by a discrete index λ, and we denote the IURs as
U λ (g).
As an example, consider the 2D rotation group SO(2).
Since SO(2) is commutative, its representation matrices
can be jointly diagonalized and so the IURs of SO(2) are
1 × 1 matrices:λU00(g) = eiλg .
(3)They satisfy the composition rule eiλ(g+h) = eiλg eiλh ,which is the manifestation of eq. 2 for this representation.
The standard Fourier series of a function on the circle is an
expansion in terms of these matrix elements, which shows1

Technically, we have defined the slightly easier to understandnotion of indecomposability, which in this context implies irreducibility.

λ

mn

√

λ
λ
(g) = dim λ Umn(g) are the L2 -normalized
where Tmn
matrix elements and η are the Fourier coefficients of f .

Integrating eq. 5 against a matrix element and using orthonormality, we find:Zλλ (g)dµ(g).
ηmn=
f (g)Tmn(6)
GThis is the (generalized) Fourier transform for compact
groups, denoted
η = Ff.
Fast and exact algorithms for the computation of Fourier
coefficients from samples of bandlimited functions on therotation groups SO(2) and SO(3) have been developed,
and the theory required to construct such algorithms forgeneral compact Lie groups is understood (Maslen &Rockmore, 1997). The group SO(2) is isomorphic to thecircle, so for G = SO(2) equation 5 reduces to a standard Fourier series on the circle, for which the well-known
O(L log L) FFT algorithm can be used. The SO(3) FFT
has complexity O(L3 log2 L) for bandlimit L (Maslen &Rockmore, 1997; Kostelec & Rockmore, 2008; Potts et al.,Harmonic Exponential Families2009). This is a tremendous speedup compared to the naive
O(L6 ) algorithm, and the algorithms presented in this paper would certainly not be feasible without the FFT.
In section 4 we discuss how these generalized FFT algorithms can be used to efficiently compute moments, but
first we discuss the generalization of the Fourier transformon compact Lie groups to the Fourier transform on certain
manifolds that are not groups.

spherical harmonics Ymλ , which form a basis for L2 (S 2 ).
Fast spherical Fourier transform algorithms were developed by Driscoll & Healy (1994).
3.4. Exponential families
An exponential family is a class of densities of the form:p(g | η) =1
exp (η · T (g)).
Zη

(8)

3.3. Harmonic analysis on homogeneous spaces
A homogeneous space for a Lie group G is a manifold Hsuch that for any two points x, y ∈ H we can find a transformation g ∈ G with gx = y. For example, the plane
is a homogeneous space for the translation group, and thesphere is a homogeneous space for the 3D rotation group.
The plane is not a homogeneous space for the 2D rotationgroup, because points at different radii cannot be rotatedinto each other.
If we pick an origin o ∈ H, such as the north pole of thesphere, we can identify any other point h ∈ H by specifying how to transform the origin to get there: h = go.
This identification will not be unique, though, if there isa nontrivial subgroup K of G containing transformations
that leave the origin invariant: K = {k ∈ G | ko = o}.
This is because if h = go, then also h = gko, so bothg and gk identify h. On the sphere for example, we cantransform the north-pole into a point h by first doing an
arbitrary rotation around the north-pole axis (which leaves
the north-pole unchanged) and then rotating the result to h.
Hence, one can think of the points h = go in a homogeneous space H as sets gK = {gk | k ∈ K} (called cosets)of group elements that are equivalent with respect to theireffect on an arbitrarily chosen origin o of H. It follows
that one can think of functions on a homogeneous space asfunctions on the group, with the special property that theyare (right) invariant to transformations from K:
f (gk) = f (g) ∀g ∈ G, k ∈ K,

(7)

because right-multiplication by k will only shuffle the elements within each coset. Finally, one can show (see supplementary material) that in a suitable basis, a subset of thematrix elements of IURs form a basis for the linear space ofsquare-integrable functions on G with this invariance property, which allows us to define the Fourier transform also
for functions on H.
The exact same equations (6 and 5) that define the Fourier
and inverse Fourier transform for a compact Lie group, define these transforms for a compact homogeneous space,
λbut only a subset of the coefficients ηmnwill be non-zero.
λFor the sphere, the matrix elements Tm0
are equal2 to the2

Various normalization and phase conventions are in use for

It is determined by a choice of sufficient statistics T , thattake the random variable g and produce a vector of real
statistics T (g).
To learn the parameters η, one can perform gradient-based
optimization of the log-likelihood of a set of iid samples
g1 , . . . , gN . The gradient is the moment discrepancy:
N

∇η

1 X
ln p(gi | η)
N i=1!
= T̄ − Ep(g|η) [T (g)],(9)

PNwhere T̄ = N1 i=1 T (gi ) are the empirical moments.
There is generally no closed form for the analytical moments (the expectation in eq. 9), so a numerical approximation is needed.

4. Harmonic Exponential Families
We define a harmonic exponential family on a group or homogeneous space as an exponential family where the sufficient statistics are given by a finite number of matrix elements of IURs. This makes sense only if the function
η · T (g) is real-valued, so that it can be interpreted as an
unnormalized log-probability. The easiest way to guarantee this is to take η to be real, and to use real functions
obtained as a sparse linear combination of complex matrix elements T (g) as sufficient statistics. For example,
1 iλθ
+ e−iλθ ) = cos(λθ). From here on, we take T
2 (eto be real, L2 -normalized functions and F the expansionof a real function in terms of real basis functions T .
The key observation that leads to an efficient algorithm forcomputing the moments of a harmonic density is that the
moments of such a density are its Fourier coefficients:
Z
Ep(g|η) [T (g)] =p(g|η) T (g) dµ(g) = F p. (10)G

Hence, one can obtain all J moments at once by sampling
p on a finite grid and then computing its Fourier transformusing a fast algorithm. As explained in section 4.1, thediscretization error can be made extremely small using only
O(J) spatial samples.
λλthe spherical harmonics, but it is enough to know that Ym∝ Tm0.

Harmonic Exponential Families

However, even evaluating p at a single position takes O(J)computations when using J sufficient statistics so that theoverall complexity is still O(J 2 ). Furthermore, in order toevaluate p we need to know the normalizing constant Zη .
The following derivation shows that we can work with theunnormalized density ϕ(g | η) = exp (η · T (g)) instead:Zλλ[Fp]mn =p(g | η) Tmn(g) dµ(g)GZ
1λϕ(g | η) Tmn(g) dµ(g)
=
(11)
Zη G
[Fϕ]λmn,=[Fϕ]000
0
The last step uses the fact that T00(g) = 1 so that [Fϕ]000
is equal to the normalizing constant:Z0
0[Fϕ]00 =ϕ(g | η) T00(g) dµ(g) = Zη(12)

the spectrum of a smooth function decays to zero asymptotically faster than O(1/λn ) for any n. So our function willbe “effectively bandlimited”, in the sense that coefficientsfor λ greater than some pseudo-bandlimit will have negligible values. If L is the maximum degree of the sufficient
statistics (the bandlimit of η · T (g)), one can obtain nearexact moments by computing the FFT up to the pseudobandlimit αL for some oversampling factor α. In practice,
we use values for α ranging from 2 to 5.

5. Harmonic Densities as Conjugate Priors
In this section we discuss the Bayesian transformation inference problem, where the goal is to infer a posterior overa Lie group of transformations given only a set of correspondence pairs (such as images before and after a camera
motion). It turns out that the harmonic exponential families are the conjugate priors for this problem, and again, thegeneralized FFT is key to performing efficient inference.

G

Next, observe that we can evaluate ln ϕ efficiently at O(J)spatial points using the inverse FFT:
ln ϕ(g | η) = η · T (g) = F −1 η (g)(13)
This computation is exact, because the log-density is bandlimited (i.e. there are only finitely many parameters).
Element-wise exponentiation then gives us ϕ evaluated ona grid.
So we have an efficient algorithm for computing moments:
1. Compute ϕ = exp (F −1 η).
2. Compute M = F ϕ
0
3. Compute Ep(g|η) [T (g)] = M/M00.

To make this computation numerically stable for highly
peaked densities, one should apply the “log-Fourier-exp”
trick described in the supplementary material.

The observed data in the Bayesian transformation inference
problem are pairs of vectors (x, y) in RD that could represent images, space-time blocks of video, point-clouds,
optical-flow fields, fitted geometric primitives, parametersof a function or other objects. In order to infer anything
about a latent transformation g, we must know the grouprepresentation R(g) that acts on the observed data. If ourdata is an image x : R2 → R, we get a representation on
the Hilbert space in which x lives: [R(g)x](p) = x(g −1 p),
where p is a point in the plane. In the finite-dimensional
analogue, where x is a vector of pixel intensities, R(g) will
be close to a permutation matrix that takes each pixel to its
proper new position. For compact groups3 this representation is unitary, and this is what we will assume for R(g)from now on. If the representation is not known in advance,
it can also be learned from data (Cohen & Welling, 2014).
If we assume that observation x is the g-transformed version of y with some independent Gaussian noise with variance σ 2 added, the likelihood function is given byp(x | y, g) = N (x | R(g)y, σ 2 ).

(14)

4.1. Approximation quality
Even though the Fourier coefficients are defined as definite integrals (eq. 6), the discrete FFT algorithms compute exact Fourier coefficients, provided the function fromwhich the discrete samples were gathered is bandlimited.
λA function is bandlimited if the coefficients ηmn
are zerofor λ greater than the band-limit L. Although the functionln ϕ(g) = η · T (g) is bandlimited, the function ϕ(g) =exp (η · T (g)) is not, so the computed coefficients are notexactly equal to the Fourier coefficients of ϕ.
However, the function ϕ(g) is smooth (infinitely differentiable), and a standard result in Fourier analysis shows thatAs discussed in section 3.2, we can bring R in blockdiagonal form by a unitary change of basis: U (g) =F −1 R(g)F . The matrix U is block diagonal, and theblocks U λ are equal to the L2 -normalized√ sufficient statistics T λ up to a scale factor: T λ (g) = dim λ U λ (g). To
simplify the computations, we shall work with data in thisnew basis: x̂ = F x so that p(x̂ | ŷ, g) = N (x̂ | U (g)ŷ, σ 2 ).
If we now choose as prior p(g) a member of the harmonic
exponential family on G, the posterior p(g | x̂, ŷ) is of the3This more generally true for unimodular groups.

Harmonic Exponential Families

same form as the prior:p(g | x̂, ŷ) ∝ p(x̂ | ŷ, g)p(g)
12∝ exp − 2 kx̂ − U (g)ŷk + η · T (g)
2σ


1 T
∝ exp
x̂ U (g)ŷ + η · T (g)
σ2
!

X
x̂λ ŷλT
λ
λ
= exp
η + √
· T (g)σ 2 dim λλ(15)i.e. we have a conjugate prior. The derivation relies on theunitarity of the representation: in expanding kx̂−U (g)ŷk2 ,we find a term kU (g)ŷk2 which is equal to kŷk2 , making
the dependence on U (g) linear (as it is in the prior).

Figure 1. Posterior distributions (top) for two correspondence
pairs (bottom).

5.1. Example: Bayesian analysis of camera rotation
To make matters concrete, we show how to compute a
posterior over the rotation group SO(3) given two images
taken before and after a camera rotation. An image is modeled as a function x : S 2 → R on the sphere, so that acamera rotation g ∈ SO(3) acts by rotating this function
over the sphere: [R(g)x](p) = x(g −1 p). We parameterize
points p ∈ S 2 as p = (ϕ, θ) for ϕ ∈ [0, 2π] and θ ∈ [0, π].
Recall from section 3.3 that we can represent p ∈ S 2 bya coset representative gp ∈ SO(3), which we parameterize
using Euler angles as gp = (ϕ, θ, 0). The transformation
gp takes the origin of the sphere to p.
It is well known4 that in this context the matrix F —defined in the previous section as the matrix that blockdiagonalizes the representation R — is given by the spherical Fourier transform F, which can be computed by an
FFT. This means that if we represent x by its Fourier coefficients x̂ = Fx, the coefficients of the rotated function
R(g)x(p) = x(g −1 p) are given by U (g)x̂, where U (g) is a
block-diagonal matrix with irreducible representations U λas blocks. As derived in the previous section, the parameters of the posterior can easily be obtained in this basis bya block-wise outer product.
Figure 1 shows the posterior p(g, | x, y) for two synthetic
spherical images x1 and x2 , and their rotations y 1 = x1
(no rotation) and y 2 = U (0, π/3, π/2)x2 . The posterior isplotted as a 3D isocontour in the ZYZ-Euler angle parameter space (α, β, γ) ∈ [0, 2π] × [0, π] × [0, 2π]. Although
this plot looks like a box, the actual manifold has the topology of a projective 3-sphere. By construction, the density
is spread through the parameter space in a way that is consistent with this topology: no discontinuities arise where
wraparound in the parameter space occurs. This is desirable, because the wraparound is not an intrinsic property ofthe manifold.
4

The derivation can be found in the supplementary material.

Due to the symmetry of these figures certain rotations cannot be distinguished from the “true” rotation, so that themodes of the posterior distributions are supported on an
entire subgroup of SO(3). Real images will not have sucha high degree of symmetry, but it will nevertheless often be
the case that a unique optimal transformation does not exist (Ma et al., 1999). Indeed, current keypoint based transformation estimation methods can easily get confused byrepeating structures in an image, such as several identical
windows on a building. Although more experimental work
is needed, our method has the theoretical advantage that besides keypoints (which form a group representation), it canmake use of parts of the image that do not allow for keypoints to be reliably placed (such as edges), while always
providing a truthful impression of the degree to which a
unique transformation or subgroup can be identified form
the data.
5.2. MAP inference
As Bayesians we are done here, but for some applications
one may wish to find a single most likely transformation.
To find the optimal transformation, first perform posterior
inference (steps 1 and 2) and then maximize (step 3):
1. Compute x̂ = F x and ŷ = F y.
1
2. Compute η̄ λ = η λ + σ2 √dimx̂ ŷ Tλ λ λ

3. Compute g ∗ = arg maxi [F −1 η̄](gi )The arg max ranges over all the points in a finite grid on
G used by the FFT synthesis F −1 . Optionally, one can refine the optimum g ∗ by performing a few steps of gradientbased optimization on η̄ · T (g) to get sub-pixel accuracy.

Harmonic Exponential Families

6. Experiments
6.1. Modelling the spatial distribution of earthquakes
We compare our model and MLE algorithm to a Kent Mixture Model (KMM) on the problem of modelling the spatial
distribution of significant earthquakes on the surface of theearth.
We obtained the Significant Earthquake Dataset (NGDC,
2015) from the National Geophysical Datacenter of the National Oceanographic and Athmospheric Administration.
In total, the dataset contains 5780 earthquakes with complete information on the position of their epicenter, and 53
earthquakes whose coordinates are missing (these were discarded in our experiments). We did not model the severityof the earthquake, but only the occurrence of significant
earthquakes (as defined by (NGDC, 2015)).

For regularization we use a diagonal Gaussian prior on η,
λwhere the precision βmcorresponding to the coefficient ofλλYm is given by βm = α dim λ = α(2λ + 1) (for someregularization parameter α). This scheme is inspired by
the fact that dim λ is the discrete Plancherel measure (Sugiura, 1990), and the empirical observation that the fitted coefficients√ become approximately uniform when weighted
λ
dim λ. Note that adding regularization does notas ηmchange the convexity of the objective function.
To find maximum a posteriori parameters η̂ for the spherical harmonic density, we perform iterative gradient-based
optimization on the log-posterior. The gradients (moment
discrepancies) are computed using the FFT-based method
described in section 4. We use the spherical FFT algorithm
implemented in the NFFT library (Keiner et al., 2009; Kunis & Potts, 2003). The gradients are fed to a standard implementation of the L-BFGS algorithm.

6.1.1. M IXTURE OF K ENT DISTRIBUTIONS
The 5-parameter Kent distribution (Mardia & Jupp, 1999)
is the spherical analogue of the normal distribution with unconstrained covariance. Being unimodal, the Kent distribution is not flexible enough to describe complicated distributions such as the spatial distribution of earthquakes. The
most flexible distribution on the sphere that we have foundin the literature is the Kent Mixture Model, first describedby Peel et al. (2001). The KMM is trained using the EM
algorithm. We use the open source Python implementationof the EM algorithm for KMMs by Höfer (2014).

6.1.3. R ESULTS
Figure 2 shows the average train and test log-likelihood
over 5 cross-validation folds, for the spherical harmonic
density and the mixture of Kent distribution. The plottedvalues correspond to the regularization settings that yielded
the best test log-likelihood.

Unlike the harmonic densities, the log-likelihood of thismodel is not convex and contains many singularities wherea mixture component concentrates on a single data pointand decreases its variance indefinitely. For this reason,
we perform randomly initialized restarts until the algorithm
has found 10 non-degenerate solutions, of which we retain
the one with the best cross-validation log-likelihood. No
regularization was used, because for the models that could
be trained within a reasonable amount of time, no overfitting was observed.

The KMM reached an average test log-likelihood of −0.37(with standard deviation of 0.03 over 5 cross-validation
folds) using 70 mixture components (5 × 70 + 69 = 419
parameters). The harmonic density reached an average test
log-likelihood of +0.3 (with standard deviation of 0.036
over 5 cross-validation folds), using bandlimit L = 140(19880 parameters). The HD still outperforms the KMM
when given a similar number of parameters: for L = 20
(440 ≈ 419 parameters), the log-likelihood is −0.28 >−0.37, with standard deviation 0.028, and for L = 19 (399parameters), the log-likelihood is −0.3 with standard deviation 0.03. The dataset and the learned densities are plottedin figure 3, clearly showing the superiority of the harmonic
density.

6.1.2. T HE S 2 HARMONIC DENSITY

7. Discussion and Future Work

The harmonic density on the 2-sphere uses spherical harmonics as sufficient statistics. The empirical moments areeasily computed using standard spherical harmonic routines, but we found that for high orders the SciPy routines are slow and numerically unstable. The supplementary material describes a simple, fast, and stable method
for the evaluation of spherical harmonics. The computation of spherical harmonics up to band-limit L = 200 (fora total of (L + 1)2 = 40401 spherical harmonics) for 5780
points on the sphere took half a minute using this method
and is performed only once for a given dataset.

What could explain the difference in log-likelihood between our model and the KMM? We believe two interrelated factors are driving this difference: the expressiveness of the model and the ease of optimization. Leaving
technicalities aside, it is clear that both the spherical harmonic exponential family and the KMM can approximate
any well-behaved density, given enough parameters and an
optimization oracle. However, as is clear from figure 2,
training time becomes prohibitive for the KMM for morethan 70 mixture components / 419 parameters, while theharmonic density can efficiently be fit for tens of thousands

Harmonic Exponential Families(a) Significant Earthquake Dataset(b) Harmonic Density

Figure 2. Top: cross-validation train (dashed line) and test (solidline) log-likelihood for the Harmonic Density (HD) and Kent
Mixture Model (KMM). Bottom: number of parameters versus
training time for both models.

of parameters.
Furthermore, the KMM training algorithm (EM) can easily
get stuck in local optima, or converge on a degenerate solution. This is the main reason for the poor runtime performance; while the KMM code could be further optimized, itis the fact that so many restarts are required to find a good
fit that makes the algorithm slow. The log-likelihood function of the harmonic density, on the other hand, is convex,
and the L-BFGS optimizer will typically converge to theglobal optimum in some 20 − 100 iterations.
As can be seen in figure 3, the harmonic density producesslight ringing artifacts that can be seen only in a log-plot
such as this. These are the result of the limited bandwidthof the log-density, and will become progressively less pronounced as the number of parameters is increased. Whilethey are clearly visible in log-space, the actual difference
between peaks and valleys is on the order of 10−3 for bandlimit L = 100. The artifacts are not visible on a nonlogarithmic plot (and in such a plot the KMM density ishardly visible at all when plotted on the same intensity
scale as the harmonic density, because the peaks are much
lower). The harmonic density also tends to prefer heavier
tails, which is probably accurate for many problems.
An interesting direction for future work is the extension
to non-compact groups. While the mathematical theory
becomes much more technical for such groups, (Kyatkin& Chirikjian, 2000) have already succeeded in developing
FFT algorithms for the Euclidean motion group which is
non-compact. From there it should be relatively straight-(c) Kent Mixture Model

Figure 3. Log-probability for harmonic density and Kent mixture
model, plotted using a perceptually accurate linear-lightness colormap on the same intensity scale.

forward to develop harmonic densities on the Euclideangroup. Harmonic densities on the Euclidean, affine or evenprojective group should find many applications in robotics
and computer vision (see e.g. (Kyatkin & Chirikjian,1999)).

8. Conclusion
We have studied a class of exponential families on compact
Lie groups and homogeneous manifolds, which we call harmonic exponential families. We have shown that for thesefamilies, maximum likelihood inference, posterior inference and mode seeking can be implemented using very efficient generalized Fast Fourier Transform algorithms. In theBayesian setting, we have shown that harmonic exponential
families appear naturally as conjugate priors in the generictransformation inference problem. Our experiments showthat training harmonic densities is fast even for very large
numbers of parameters, and that far superior likelihood canbe achieved using these models.

Acknowledgements
This research was supported by NWO (grant number
NAI.14.108), Facebook and Google.

References
Beran, Rudolf. Exponential Models for Directional Data.
The Annals of Statistics, 1979.
Chirikjian, G.S. and Kyatkin, A.B. Engineering Applications of Noncommutative Harmonic Analysis. CRC

Harmonic Exponential Families

Press, 1 edition, May 2001. ISBN 9781420041767.
Cohen, T. and Welling, M. Learning the Irreducible Representations of Commutative Lie Groups. In International
Conference on Machine Learning (ICML), volume 32,
2014.
Diaconis, P. Group representations in probability andstatistics. Hayward, CA: Institute of Mathematical
Statistics, 1988. ISBN 0940600145.
Driscoll, JR and Healy, DM. Computing Fourier transforms and convolutions on the 2-sphere. Advances inapplied mathematics, 1994.
Gatto, Riccardo and Jammalamadaka, Sreenivasa Rao. The
generalized von Mises distribution. Statistical Methodology, 4(3):341–353, 2007. ISSN 15723127.
Goodman, Roe and Wallach, Nolan R. Symmetry, Representations, and Invariants. Springer, 2009. ISBN 9780-387-79852-3.
Höfer, Sebastian. Kent mixture model for scikit-learn,
2014. URL https://github.com/shoefer/scikit-learn.
Keiner, J., Kunis, S., and Potts, D. Using NFFT 3 - asoftware library for various nonequispaced fast Fouriertransforms. ACM Transactions on Mathematical Software, 36:1–30, 2009.
Kondor, Risi. Group theoretical methods in machine learning. PhD thesis, Columbia University, 2008.
Kostelec, Peter J. and Rockmore, Daniel N. FFTs on therotation group. Journal of Fourier Analysis and Applications, 14(2):145–179, 2008.
Kunis, Stefan and Potts, Daniel. Fast spherical Fourier algorithms. Journal of Computational and Applied Mathematics, 161:75–98, 2003.
Kyatkin, AB and Chirikjian, GS. Computation of Robot
Configuration and Workspaces via the Fourier Transform on the Discrete-Motion Group. The International
Journal of Robotics Research, 1999.
Kyatkin, Alexander B and Chirikjian, Gregory S. Algorithms for Fast Convolutions on Motion Groups. Applied
and Computational Harmonic Analysis, 9(2):220–241,
September 2000.
Ma, Y, Soatto, S, Kosecka, J, Sastry, S, and Košecká, J. Euclidean reconstruction and reprojection up to subgroups.
In The Proceedings of the Seventh IEEE International
Conference on Computer Vision (ICCV), volume 2, pp.
773—-780. California Univ., Berkeley, CA, IEEE, 1999.

Mardia, K. V. and Jupp, P. E. Directional Statistics. John
Wiley & Sons, 1 edition, 1999. ISBN 0471953334.
Maslen, DK and Rockmore, DN. Generalized FFTs – a
survey of some recent results. Groups and Computation
II, 1997.
NGDC. National Geophysical Data Center / World Data
Service: Significant Earthquake Database. [accessed 0118-2015], 2015. URL http://www.ngdc.noaa.
gov/nndc/struts/form?t=101650&s=1&d=1.
Peel, David, Whiten, William J, and McLachlan, Geoffrey J. Fitting Mixtures of Kent Distributions to Aid in
Joint Set Identification. Journal of the American Statistical Association, 96(453):56–63, 2001. ISSN 0162-1459.
doi: 10.1198/016214501750332974.
Pinchon, Didier and Hoggan, Philip E. Rotation matrices
for real spherical harmonics: general rotations of atomic
orbitals in space-fixed axes. Journal of Physics A: Mathematical and Theoretical, 40(7):1597–1610, February
2007. ISSN 1751-8113.
Potts, Daniel, Prestin, J, and Vollrath, A. A fast algorithm for nonequispaced Fourier transforms on the rotation group. Numerical Algorithms, pp. 1–28, 2009.
Sugiura, Mitsuo. Unitary Representations and Harmonic
Analysis. John Wiley & Sons, New York, London, Sydney, Toronto, 2nd edition, 1990.
Wainwright, Martin J. and Jordan, Michael I. Graphical Models, Exponential Families, and Variational Inference. Foundations and Trends R in Machine Learning,
1(1–2):1–305, 2007. ISSN 1935-8237.