Directional Statistics with the Spherical Normal
Distribution
Søren Hauberg
Department of Mathematics and Computer Science
Technical University of Denmark
Kgs. Lyngby, Denmark
sohau@dtu.dk

Abstract—A well-known problem in directional statistics — thestudy of data distributed on the unit sphere — is that current
models disregard the curvature of the underlying sample space.
This ensures computationally efficiency, but can influence results.
To investigate this, we develop efficient inference techniques fordata distributed by the curvature-aware spherical normal distribution. We derive closed-form expressions for the normalization
constant when the distribution is isotropic, and a fast and accurate approximation for the anisotropic case on the two-sphere. We
further develop approximate posterior inference techniques forthe mean and concentration of the distribution, and propose a fastsampling algorithm for simulating the distribution. Combined,
this provides the tools needed for practical inference on the unit
sphere in a manner that respects the curvature of the underlying
sample space.

I. I NTRODUCTION
Directional statistics considers data distributed on the unit
sphere, i.e. kxn k = 1. The corner-stone model is the vonMises-Fisher distribution [19] which is derived by restricting
an isotropic normal distribution to lie only on the unit sphere: κ

vMF(x | µ, κ) ∝ exp − kx − µk2 ∝ exp (κ x| µ) ,2
kxk = kµk = 1, κ > 0.
(1)
Evidently, the von Mises-Fisher distribution is constructed fromthe Euclidean distance measure kx − µk, i.e. the length of theconnecting straight line. A seemingly more natural choice ofdistance measure is the arc-length of the shortest connecting
curve on the sphere, which amounts to the angle between two
points, i.e. dist(x, µ) = arccos(x| µ). From this measure, it
is natural to consider the distribution


λ
2 |
SN(x | µ, λ) ∝ exp − arccos (x µ) ,2(2)
kxk = kµk = 1, λ > 0,with concentration parameter λ. This is an instance of the
Riemannian normal distribution [20] which is naturally deemed
the spherical normal distribution. This is the topic of interestin the present manuscript.
At this point, the reader may wonder if the choice ofdistance measure is purely academic or if it influences practical
This project has received funding from the European Research Council(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement n° 757360). SH was supported by a research
grant (15334) from VILLUM FONDEN.

Fig. 1. Likelihood of the spherical mean under the von Mises-Fisher andspherical normal models. Top row: Two observations are moved further andfurther apart until they are on opposite poles (coloring show the sum of thelikelihood terms for the data). Middle row: The likelihood of the mean under avon Mises-Fisher distribution. As the data moves further apart, the variance ofthe likelihood grows until it degenerates into the uniform distribution. Bottom
row: The likelihood of the mean under a spherical normal distribution. As
the observations moves apart, the likelihood becomes increasingly anisotropic
until it stretches along the equator. This is the most natural result.

inference. We therefore consider a simple numerical experiment,
where we evaluate the likelihood of the mean µ of a von MisesFisher and a spherical normal distribution. We consider one
observation on the north pole, and another observations which
we move along a great circle from the north to the south pole.
When the two observations are on opposite poles the mean µ
cannot be expected to be unique; rather intuition dictate that any
point along the equator (the set of points that are equidistantto the poles) is a suitable candidate mean. Figure 1 shows theexperiment along with the numerically evaluated likelihood ofµ. We observe that under the von Mises-Fisher distribution,
this likelihood is isotropic with an increasing variance as theobservations move further apart. When the observations areon opposite poles, the likelihood becomes uniform over theentire sphere, implying that any choice of mean is as good as
another. This is a rather surprising result, that align poorly withgeometric intuition. Under the spherical normal distribution,
the likelihood is seen to be anisotropic, where the varianceincrease more orthogonally to the followed great circle, thanit does along the great circle. Finally, when the observation
reaches the south pole, the likelihood concentrates in a “belt”

Fig. 2. Flat versus curved metrics. Straight lines (purple) correspond to
Euclidean distances used by the von Mises-Fisher distribution. This flat metricis in contrast to the spherical arc-length distance (yellow geodesic curves).
Under this distance measure, geodesic triangles are “fatter” than corresponding
Euclidean triangles. This distortion, which is the key characteristics of thesphere, is disregarded by the von Mises-Fisher.

Fig. 3. The tangent space Tµ of the sphere and related operators.

logarithm map,
Logµ (x) = (x − µ (x| µ))|θ
,sin(θ)(3)

θ = arccos(x µ),
along the equator, implying that any mean along the equator is0
a good candidate. This likelihood coincide with the geometric with the convention /sin(0) = 1. The inverse mapping, known
as the exponential map, moves a tangent vector v back to theintuition.
sphere,
The core issue is that the Euclidean distance measure appliedby the von Mises-Fisher distribution is a flat metric [8] as
Expµ (v) = µ cos(kvk) + sin(kvk)/kvkv.
(4)
illustrated in Fig. 2. This imply that the von Mises-Fisher
distribution is incapable of reflecting the curvature of the sphere. A tangent vector v ∈ Tµ can be moved to another tangent
This amounts to applying a
As curvature is one of the defining mathematical properties space Tµ0 by a parallel transport.
0
rotation
R
that
move
µ
to
µ
[9].
of the spherical sample space, it is a serious limitation of theb) Distributions on the sphere.: Since the unit sphere is
von Mises-Fisher model that this is disregarded as it may giveacompact
space, we can define a uniform distribution withmisleading results. For practical inference, there is, however,densityno viable alternative that respect the curvature of the sphericalΓ(D/2)
sample space. We fill this gap in the literature and contribute
Uniform(x) = A−1
.
(5)
D−1 =
2π D/2• closed-form expressions for the normalization constant of
the isotropic spherical normal distribution (Sec. III-A) and The corner-stone distribution on the unit sphere is the alreadyan efficient and accurate approximation for anisotropic discussed von Mises-Fisher distribution [19], which is derived
by restricting the isotropic normal distribution to the unit
case on the two-sphere;
sphere. Consequently it is defined with respect to the Euclidean• maximum likelihood (Sec. IV) and approximate Bayesian
distance measure. The mean parameter of the von Misesinference techniques (Sec. V) for the spherical normal;
Fisher distribution can be estimated with maximum likelihood• an efficient sampling algorithm for the spherical normal
in closed-form, but the concentration parameter κ must be(Sec. VI).
estimated using numerical optimization [19]. Bayesian analysis
As examples of these techniques, we provide a new clustering is simplified since the distribution is conjugate with itself for themodel on the sphere (Sec. IV-C) and a new algorithm for mean, and with the Gamma distribution for the concentration.
Kalman filtration on the sphere (Sec. V-B). Relevant source
The von Mises-Fisher distribution can be extended to becode will be published alongside this manuscript.
anisotropic by restricting an anisotropic normal distribution
to the unit sphere [16], giving the Fisher-Bingham (or Kent)
distribution,
II. BACKGROUND AND RELATED WORK1p(x) = exp κγ |1 x + β[(γ |2 x)2 − (γ |3 x)2 ] ,(6)Ca) Spherical geometry.: We start the exposition with a where the normalization constant C can be evaluated in closedbrief review of the geometry of the unit sphere. This is a form on S 2 . Like the von Mises-Fisher distribution, this use
compact Riemannian manifold with constant unit curvature. It the Euclidean distance function, and therefore does not respect
2π D/2has surface area AD−1 = Γ(the curvature of the sphere.
D/2) , where Γ is the usual Gamma
function. It is often convenient to linearize the sphere around a
One approach to constructing an anisotropic curvaturebase point µ, such that points on the sphere are represented in aware distribution is to assume that Logµ (x) follow a zeroa local Euclidean tangent space Tµ (Fig. 3). The convenience mean anisotropic normal distribution [28]. This tangential
stems from the linearity of tangent space, but this comes at normal distribution is simple, and since it is defined over thethe cost that Tµ gives a distorted view of the curved sphere. Euclidean tangent space, standard Bayesian analysis hold forSince the maximal distance from µ to any point is π, we are its precision (i.e. it is conjugate with the Wishart distribution).
generally only concerned with tangent vectors v ∈ Tµ where However, the distribution disregard the distortion of the tangent
kvk ≤ π. A point x ∈ S D−1 can be mapped to Tµ via the space: following the change of variable theorem we seethat the tangential normal over S 2 has density p(x) =N (x | µ, Σ) |det(∂x Logµ )| = N (x | µ, Σ) θ/sin(θ), whereθ = arccos(x| µ). Surprisingly, this distribution is thus bimodal
with modes at ±µ. This rather unintuitive property indicatethat one should be careful when applying this model.
An approach related to ours is that of Purkayastha [21]who study p(x) ∝ exp(−κ arccos(x| µ)) and show that the
maximum likelihood estimate of µ is the intrinsic median on
the sphere.
c) Manifold statistics.: The spherical normal distribution,
which we consider in this manuscript, is an instance of thegeneral Riemannian normal distribution. While this distribution
has seen significant theoretical studies [20, 29], practical tools
for inference are lacking. Generally, its normalization constantdepend on both mean and concentration parameters and isunattainable in closed-form. Even evaluating the gradient of
the log-likelihood with respect to the mean, thus, require
expensive Monte Carlo schemes [5, 31]. Similar concerns holdfor the concentration parameters. A key contribution of thispaper is that we provide tools for the spherical setting thatdrastically simplify and speed up parameter estimation. We
further provide tools for approximate Bayesian inference on
Riemannian manifolds — a topic that has yet to be addressedin the literature.
d) Applications of directional statistics.: The use of
directional distributions span many branches of science, ranging
from modeling wind directions [19], protein shapes [10], gene
expressions [6] to fMRI time series [22]. In these applications,the directional model appear due to normalization of thedata. Histograms, such as word counts, can also be efficiently
modeled as directional data [6]; in fact, the arc-length distanceon the sphere correspond to the distance between discrete
distributions under the Fisher-Rao metric [4], which is anotherindication that we need curvature-aware models. Other recent
applications include clustering surface normals to aide robot
navigation [27, 28], and speaker direction tracking [30]. At
the foundational level, directional statistics also appear inprobabilistic estimates of the principal components of Euclidean
data [14, 26].
III. BASIC PROPERTIES OF THE SPHERICAL NORMAL
The spherical normal distribution as presented in Eq. 2 isisotropic. Noting that arccos2 (x| µ) = Logµ (x)| Logµ (x) is
merely the squared Euclidean norm measured in the tangent
space Tµ allows us to generalize Eq. 2 to be anisotropic:
1
SN(x | µ, Λ) ∝ exp − Logµ (x)| ΛLogµ (x) ,(7)2where the concentration parameter Λ is a symmetric positive
semidefinite matrix. This is defined in Tµ and consequently
has an eigenvector µ with eigenvalue zero.
When the concentration goes towards zero, it is easy to seethat the spherical normal tends towards the uniform distribution,i.e. when Λ = λI we havelim SN(x | µ, λI) = Uniform(x).

λ→0+

(8)

When λ goes towards infinity, the spherical normal becomes adelta function. Similar properties hold for the von Mises-Fisher
distribution [19].
It is worth noting that this general spherical normal distribution is the spherical distribution with maximal entropy for agiven mean and covariance [20]. Note that this is with respect
to the probability measure induced by the spherical metric.
Similar statements hold for the von Mises-Fisher distribution[19], but with respect to the measure induced by the Euclidean
metric.
The spherical normal, thus, share many desirable properties
with the von Mises-Fisher distribution, while still being
curvature-aware.
A. Normalizing the spherical normal
Equation 2 only specifies the spherical normal distribution
up to a constant. In the online supplements1 we show thefollowing result:
Proposition 1 (Isotropic normalization). The normalization
constant of the isotropic spherical normal distribution (2) over
S D−1 is given by Eq. 9a when D is even and by Eq. 9b when
D is odd (see Fig. 4). Here erf is the imaginary error function[1], and Re[·] and Im[·] takes the real and imaginary parts ofa complex number, respectively.
Remark 1. While the complex unit appear in Eq. 9 the entire
expression evaluates to a real number.
The anisotropic spherical normal distribution (7) does,
unfortunately, not appear to have a closed-form expression
for its normalization constant, and we must resort to approximations. Here we consider the anisotropic spherical normal
over the ever-present S 2 and provide a simple deterministic
approximation that, in our experience, provides a relative errorof approximately 10−4 . This is more than sufficient for practicaltasks. To derive this, we first express the distribution in thetangent space of the mean,Z
1Z2 (Λ) =exp − Logµ (x)| ΛLogµ (x) dx
(10)
2
2


ZS
1
=
exp − v| Λv det(J)dv,
(11)
2
kvk≤π
where v = Logµ (x) and det(J) = sin(kvk)/kvk denotes theJacobian of Expµ (x). This expression is independent of the
choice of orthogonal basis of the tangent space, so we need
only consider the case where Λ = diag(λ1 , λ2 , 0). WritingEq. 11 in polar coordinates (r, θ) then shows
Z 2π Z πn r2
Z 2 (Λ) =exp − (λ1 cos2 (θ)
2
0
0
(12)
o
2
+ λ2 sin (θ)) sin(r)drdθ
Z 2π

1
=
(13)
Z1 (λ1 − λ2 ) cos2 (θ) + λ2 dθ.
2π 0
1 http://www2.compute.dtu.dk/~sohau/papers/fusion2018/hauberg_
fusion2018_supplements.pdf

r !
r
r
D
D−2
AD−2 2π (−1) /2−1
π
λ
erf
π
+
D−3
D/2 − 1
2
2λ
2
2
λ
2
D/2−2



  

X
D−2
(D − 2 − 2k)2πλ − (D − 2 − 2k)i
√
·
(−1)k
exp −
Re erf
k
2λ
2λ

AD−2
Z1(even) (λ) = D−2



(9a)

k=0

(D − 3)/2(D − 3)/2





D−2
(D − 2 − 2k)2
exp −
2
2λ
k
  
k=0  (D − 2 − 2k)iπλ − (D − 2 − 2k)i√√
· Im erf+ Im erf
2λ
2λ

(−1)
Z1(odd) (λ) = AD−2
D−3

r

π
2λ

X

(−1)k

(9b)

Fig. 4. Normalization constants for the isotropic spherical normal distribution, when D is even and odd. Here erf is the imaginary error function, while Re[·]and Im[·] takes the real and imaginary parts of a complex number, respectively.

where R is a rotation matrix. Two independent sphericalnormally distributed variables can (informally) be “added” by
convolving their distributions. When these variables have thesame mean, thenpx (x) = SN(x | µ, Λx )⇒py (y) = SN(y | µ, Λy )(17)
(px ? py )(z) = SN (z | µ, Λx + Λy ) ,Fig. 5. Left: The inverse of Z1 (λ). Right: Relative errors when computing
Z2 (Λ).

This does not appear to have a closed-form expression, butdoes lend itself to a good approximation. Noting that 1/Z1 (λ) isalmost a straight line (see left panel of Fig. 5), we approximate
Z 2π
1Z2 (Λ) ≈dθ(14)
2 (θ) + λ ) + b
a
((λ
−
λ
)
cos
1
2
2
0
1
=p
,
(15)
det (aΛ + bI)

where ? denotes convolution. In Sec. V-B we use both results
in a spherical Kalman filter.
C. A note on covariances
The dispersion of the spherical normal is parametrized by aconcentration matrix Λ. From this the variance and covarianceof the distribution is defined as [20]Z
Var[x] =(18)
arccos2 (x| µ) SN(x | µ, Λ)dx
D−1
SZ
Cov[x] =Logµ (x)Logµ (x)| SN(x | µ, Λ)dx.
S D−1

The empirical counterparts of these expressions are found by
where a and b are estimated by fitting a straight line to 1/Z1 (λ). replacing integrals with averages as usual. Generally we can
In Sec. IV-C we consider an EM algorithm for mixtures of evaluate these expressions numerically by sampling (Sec. VI),
spherical normal distributions. Here we track the relative though we note that in the case of an isotropic spherical normal
approximation error of the normalization constant (using over S 2 , the variance can be expressed in closed-form (seeexpensive numerical integration as ground truth) through the supplements1 ). The expression is, however, somewhat moreentire run of the EM algorithm. The right panel of Fig. 5 shows convoluted than in the Euclidean setting where the variance is
a histogram of these errors. From this data, we observe that the inverse concentration (precision). Note that the variance isthe relative error of our proposed approximation is between bounded by π 2 since the maximal spherical distance is π.
2.1 · 10−9 and 6.1 · 10−4 with an average of 7.9 · 10−5 . This
IV. M AXIMUM LIKELIHOOD ESTIMATION
is plenty accurate for practical inference tasks.
We now consider maximum likelihood estimation for the
spherical normal. As usual, we write the log-likelihood of data
B. Rotations and convolutions
x1:N
In the spherical domain we cannot perform addition and!
NYmultiplication, but similar operations are available. Here wef (µ, Λ) = log
SN(xn | µ, Λ)
state without proof that spherical normals are closed under(19)
n=1
rotationN
1X=−
Logµ (xn )| ΛLogµ (xn ) − N log(Z2 (Λ)).
x ∼ SN(µ, Λ) ⇒ Rx ∼ SN(µ, RΛR| ),
(16)
2
n=1

This initialization coincides with the standard “least squares”
estimator of the inverse covariance on general Riemannian
manifolds [5]. It is therefore not surprising that only few
iterations are required, and when computational speed is of
the utmost importance, this estimator is a good approximation.
Pennec [20] provides an alternative approximate estimator forsmall concentration matrices that compensates for the curvatureof the sample space.

Algorithm 1 Offline maximum likelihood: µ
1: Initialize µ as a random data point.
2: repeat
PN
3:∇µ ← − N1 n=1 Logµ (xn ).
4:µ ← Expµ − 12 ∇µ .
5: until k∇µ k < 10−5 .

Algorithm 2 Online maximum likelihood: µ

C. Example: EM algorithm for mixture models

1: Initialize µ ← x1 .
2: for n = 2 to . . . do3:µ ← Expµ n1 Logµ (xn ) .
4: end for

A. The mean
Since the normalization constant does not depend on µ, it
is easy to see that the derivative of the log-likelihood withrespect to the mean is
N
∂f
1 X
=−
Logµ (xn )| Λ.
∂µ
N n=1

(20)

The mean can then be found by standard Riemannian gradientdescent [2], but it is easier
PN to consider the steepest descent [5]given by ∇µ = − N1 n=1 Logµ (xn ). The optimization thenamounts to mapping the data to the tangent at the current µ,
computing its Euclidean average, and mapping it back to thesphere. This is summarized in Algorithm 1.
Alternatively, the mean can be computed by repeatedgeodesic interpolation [23], which can be done in an online
fashion. From two points, the mean is estimated as the midpointof the connecting geodesic; when given a third point, themean is updated by going 1/3 along the geodesic between
the two-point mean and the third point, and so forth. This issummarized in Algorithm 2. Salehian et al. [23] show that thissimple estimator converges to the true mean. This algorithm canalso be seen as stochastic gradient descent with a particularly
efficient step-size selection.
B. The concentration matrix
Using our approximation of the normalization constant (15),the gradient of the log-likelihood with respect to Λ is (expressedin a two-dimensional basis of Tµ )
( N
∂f
1 X
=−
Logµ (xn )Logµ (xn )|
∂Λ
2 n=1
)
(21)
− aN (a Λ + b I)

−1.

It is straight-forward to extend the presented maximum
PK likelihood estimators to mixture models, i.e. p(x) =Pk=1 πk SN(x | µk , Λk ), where the responsibilities πk satisfyk πk = 1, πk ≥ 0. The derivation follows the well-known
Euclidean analysis [7] and will not be repeated here.
As an example, we consider clustering of surface normals
extracted from depth images [11, 25]. This is, e.g., useful for
robot navigation [27, 28]. In practice, such data is highly noisy
due to limitations of the depth-sensing camera, so we extend the
mixture model toPinclude a component with uniform distribuKtion, i.e. p(x) = k=1 πk SN(x | µk , Λk )+πK+1·Uniform(x).
The left panel of Fig. 6 shows example data, the contours of
the estimated spherical normals, and the corresponding regions
on the sphere where each component is dominant. The center
panels show the corresponding image segmentation. This seems
to correspond well to the scene geometry. It is interesting tonote that the uniform component mostly captures edge andshadow areas where surface normals are unstable. In the figure,
we consider 5 spherical normal components.
Next we sample 50 random depth images from the NYU
Depth Dataset [25], which depicts indoor environments. In
such scenes, a few surface normals tend to dominate as mostfurniture has similar normals as the walls, floors and ceilings.
We perform clustering on these 50 images using both a mixtureof spherical normal distributions and a mixture of von MisesFisher distributions. For each model, we select the optimal
number of components using the Akaike Information Criteria(AIC) [3] measured on held-out data. The right-most panel ofFig. 6 show a histogram of the number of selected components
for both mixture models. We see that significantly fewer
components are needed when using a mixture of sphericalnormals than when using a mixture of von Mises-Fisher
distributions. This indicate that the von Mises-Fisher model
tend to over-segment the depth images.
V. A PPROXIMATE BAYESIAN INFERENCE
Thus far, we have considered maximum likelihood estimation
for the spherical normal, however, some may prefer a Bayesianmachinery. Here we consider the mean and concentration
separately.
A. Unknown mean, known concentration

An optimum can be found using standard Riemannian optimization [2] over the symmetric positive definite cone.
Initializing with the inverse covariance of Logµ (x) usually
ensures convergence in 3−4 iterations.

Setting a spherical normal prior over the mean p(µ) =SN(µ | µ0 , Λ0 ) while
Q assuming a spherical normal likelihoodp(x1:N | µ, Λ) = n SN(xn | µ, Λ) quickly reveals that theresulting posterior is not a spherical normal. This suggests

Fig. 6. Clustering of surface normals; see text for details.

that approximations are in order. First, we note that the logposterior,g(µ) = log(p(µ | x1:N , Λ, µ0 , Λ0 ))
=−

(22)

N
X

1
Logµ (xn )| ΛLogµ (xn )
2 n=1

(23)

1
− Logµ0 (µ)| Λ0 Logµ0 (µ) + const,2
is simultaneously expressed in tangent spaces at both µ andµ0 , which complicates the analysis. Due to the symmetry of
the spherical normal, we can, however, rewrite Eq. 23 asNg(µ) = −1X
Logµ (xn )| ΛLogµ (xn )
2 n=1

(24)

1
− Logµ (µ0 )| RΛ0 R| Logµ (µ0 ) + const,2
where R is the parallel transport from Tµ0 to Tµ . The maximum
a posteriori (MAP) estimate of the mean can then be foundwith Riemannian gradient descent, where the gradient is
(25)
N
∂g
1X
1
|
|
|
=−
Logµ (xn ) Λ − Logµ (µ0 ) RΛ0 R .
∂µ
2 n=1
2

yt be the unobserved variable at time t with initial density
y0 ∼ SN(µ0 , Λ0 ) and assume the predictive distribution (26)p(yt+1 |yt ) = SN yt+1 Rµ, RΛ0 R| + [Λpred ]µ ,where R is a rotation matrix encoding the expected temporaldevelopment and [Λpred ]µ denotes the concentration of thepredictive distribution expressed in the basis of Tµ . Assuming a spherically normal likelihood, the suggested Laplaceapproximation can be used to estimate the posterior.
We implement this for tracking the direction of the left femur(thigh bone) of a person walking [12]. Figure 7B shows theoriginal data (blue), the spherical normal filtered path (orange),
and a von Mises-Fisher filtered path (purple). We see that thespherical normal filter smooths the observed data as expected,
while the von Mises-Fisher filter appears to have a bias thattakes the filtered path outside the data support.
C. Unknown concentration, known mean
Following standard approaches, we here assume a Wishartprior p(Λ | Λ0 , m) = W(Λ | Λ0 , m) for the concentration,
where we assume the precision is expressed with respect toan orthogonal basis of Tµ . As before, we note that this prior
is not conjugate to the spherical normal, and again employ a
Laplace approximation of the posterior. The sample space ofthe posterior is the cone of positive definite matrices, so the
Laplace approximation need to be adapted to this constraint.
The log-posterior in this setting is1h(Λ) = m−(D−1)−1log det(Λ) − Tr Λ−1
0 Λ
2
2
(27)
N
X
1
Logµ (xn )| ΛLogµ (xn ) − N log(Z2 (Λ)).
−
2 n=1

The optimal µN can be found akin to Algorithm 1. This is,however, only a point estimate and does not approximate thefull posterior. Here we consider a Laplace approximation [7]as these generalize easily to the spherical setting. Since themetric in Tµ reduces to the identity around the origin, thisamounts to approximating the concentration of the posteriorΛN with the Riemannian Hessian of g(µ) at µN [2].
Figure 7A shows an example with a single observation witha spherical normal likelihood, as well as a spherical normal As in the maximum likelihood setting, we apply the approxiprior over the unknown mean. The intensity of the sphere is mate normalization constant (15). The derivative then become(
a numerical evaluation of the true posterior, while the orange
1∂hcurve outlines the Laplace approximation of the mean posterior.
=−(D − m)Λ−1 + Λ−10∂Λ2This appears to be a good approximation.
(28)
)
N
X
−1
B. Example: directional Kalman filter+
Logµ (xn )Logµ (xn )| − aN (a Λ + b I).
The tools presented thus far allow us to build a model akin to
the Kalman filter [15] over the unit sphere. Previous algorithms
for Riemannian Kalman filters predominantly rely on unscented
transforms to adapt to the nonlinear state space [13, 17]. Let

n=1This can then be optimized using standard Riemannian optimization [2] to find a maximum a posteriori estimate of
Λ. As for the mean, we can build a Laplace approximation

Fig. 7. A: Laplace approximation of the mean posterior. B: data (blue) and the filtered path (orange). C: samples from a Laplace approximation to theconcentration posterior. D: acceptance rates for the sampling algorithm.

Algorithm 3 Sampling x ∼ SN(µ, Λ).
1: repeat
2:
3:


−1 
I
v ∼ N 0, Λ + D−2.
πr←sin(kvk) D−2
kvk
1 |− 2 v Λ+ D−2
π I v

exp(− 12 v| Λv)(

)
.
) )

and 1/λ2 in the range (0, π 2 ]. On average the acceptance rateover this domain is 82%; for large concentrations (the common
scenario) the acceptance rate is higher. This is quite efficient.
VII. S UMMARY AND OUTLOOKThis paper contributes the first practical tools for efficient
curvature-aware inference over the unit sphere. We provide
a closed-form expression for the normalization constant ofisotropic spherical normal distributions in any dimension.
From this we build good approximations for the anisotropic
case in S 2 . We provide efficient algorithms for maximumof the posterior of Λ. Extending the standard derivation of likelihood estimation of the mean and concentration parameters
the Laplace approximation [18] shows that the approximate of the distribution, and exemplify these with an EM algorithm
posterior of Λ follow a matrix log-normal distribution over the for mixtures of spherical normals. We further provide tools
positive definite cone [24], where the concentration is given by for approximate Bayesian inference in the form of Laplace
the Riemannian Hessian of the log-posterior at ΛMAP . As an approximations for the posteriors of both the mean and theexample, Fig. 7C shows samples from the Laplace posterior concentration parameters. We exemplify these approximations
estimated from 50 observations and a Wishart prior. The choice with a spherical Kalman filter. Finally, we provide a simple
of a Wishart prior is not essential as there are no conjugate yet efficient sampling algorithm for simulating the sphericalnormal.
properties to exploit.
Several questions, however, remain open. Most importantly,
our approximation to the anisotropic normalization constant
D. Unknown mean and concentration
2
The concentration matrix Λ is fundamentally tied to the does not extend beyond S , which is a strong limitation. This
mean µ since it is expressed with respect to a basis of Tµ . is similar to established models in directional statistics, where
distribution only has a known
Changing µ, thus, renders Λ meaningless. This complicate a the anisotropic Fisher-Bingham
2
normalization
constantforS.
joint model of both µ and Λ, and we do not investigate the
The proposed Laplace approximations carry over to otherissue further in this manuscript.
Riemannian manifolds, where Bayesian inference is rarely
applied. Perhaps more importantly, the Laplace approximation
VI. S AMPLINGof the posterior concentration matrix also carry over to the
In many computational pipelines it is essential to be able
Euclidean domain, where it can be of value when approximating
to draw samples from the applied distributions. We therefore
posteriors over precision matrices. The resulting matrix logprovide a simple, yet efficient, algorithm for simulating thenormal distribution over the positive definite cone is signifispherical normal. Expressed in the tangent space of µ the
cantly more informative than the commonly applied Wishart
spherical normal has distribution
distribution, as it can capture anisotropic distributions over the


D−2
positive definite cone.
1 |sin(kvk)
SNTµ (v) ∝ exp − v Λv.
(29)
It may be beneficial to consider variational approximations to2
kvk
the posteriors rather than the proposed Laplace approximations.
D−2
We note that exp(−(D − 2)kv|k2/2π) ≥ ( sin(kvk)/kvk)
for This should be feasible when using the proposed approximate
kvk ≤ π. Since this envelope is reasonably tight, we propose normalization constant.
a basic rejection sampling [7] scheme where tangent vectors
are proposed from N (0, Λ + (D − 2)/πI). This is summarizedin Algorithm 3 and Fig. 7D shows the acceptance rate for 1/λ1
exp(

(
4:
u ∼ Uniform(0, 1).
5: until kvk ≤ π and u ≤ r.
6: x ← Expµ (v).

R EFERENCES[1] M. Abramowitz and I. A. Stegun. Handbook of mathematical functions: with formulas, graphs, and mathematical
tables, volume 55. Courier Corporation, 1964.
[2] P.-A. Absil, R. Mahony, and R. Sepulchre. Optimization
Algorithms on Matrix Manifolds. Princeton University
Press, 2008.
[3] H. Akaike. A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6):716–723, 1974.
[4] S. Amari and H. Nagaoka. Methods of informationgeometry, volume 191. American Math. Soc., 2007.
[5] G. Arvanitidis, L. K. Hansen, and S. Hauberg. A locally
adaptive normal distribution. In Advances in Neural
Information Processing Systems, 2016.
[6] A. Banerjee, I. S. Dhillon, J. Ghosh, and S. Sra. Clustering on the unit hypersphere using von mises-fisher
distributions. Journal of Machine Learning Research, 6:1345–1382, 2005.
[7] C. M. Bishop. Pattern Recognition and Machine Learning(Information Science and Statistics). Springer-Verlag New
York, Inc., Secaucus, NJ, USA, 2006.
[8] M. R. Bridson and A. Haefliger. Metric spaces of nonpositive curvature, volume 319. Springer Science &Business Media, 2011.
[9] M. do Carmo. Riemannian Geometry. Mathematics(Boston, Mass.). Birkhäuser, 1992.
[10] T. Hamelryck, J. T. Kent, and A. Krogh. Sampling realistic
protein conformations using local structural bias. PLoS
Computational Biology, 2(9):e131, 2006.
[11] M. A. Hasnat. Unsupervised 3D image clustering andextension to joint color and depth segmentation. PhD
thesis, Université Jean Monnet - Saint-Etienne, Oct. 2014.
[12] S. Hauberg. Principal Curves on Riemannian Manifolds.
IEEE Transactions on Pattern Analysis and MachineIntelligence (TPAMI), 2016.
[13] S. Hauberg, F. Lauze, and K. S. Pedersen. Unscented
kalman filtering on riemannian manifolds. Journal of
Mathematical Imaging and Vision (JMIV), 46(1):103–120,
2013.
[14] S. Hauberg, A. Feragen, R. Enficiaud, and M. J. Black.
Scalable robust principal component analysis using grassmann averages. IEEE Transactions on Pattern Analysis
and Machine Intelligence (TPAMI), 2016.
[15] R. E. Kalman. A new approach to linear filtering andprediction problems. Transactions of the ASME–Journalof Basic Engineering, 82(Series D):35–45, 1960.
[16] J. T. Kent. The Fisher-Bingham Distribution on theSphere. Journal of the Royal Statistical Society. Series B(Methodological), 44(1):71–80, 1982.
[17] G. Kurz, I. Gilitschenski, and U. D. Hanebeck. Recursive
nonlinear filtering for angular data based on circular
distributions. In American Control Conference (ACC),
2013, pages 5439–5445. IEEE, 2013.
[18] D. J. MacKay. Information theory, inference and learningalgorithms. Cambridge university press, 2003.
[19] K. V. Mardia and P. E. Jupp. Directional Statistics. Wiley,
1999.
[20] X. Pennec. Intrinsic Statistics on Riemannian Manifolds:
Basic Tools for Geometric Measurements. Journal of
Mathematical Imaging and Vision (JMIV), 25(1):127–154,
2006.
[21] S. Purkayastha. A rotationally symmetric directional
distribution: obtained through maximum likelihood characterization. Sankhyā: The Indian Journal of Statistics,
Series A, pages 70–83, 1991.
[22] R. E. Røge, K. H. Madsen, M. N. Schmidt, and M. Mørup.
Infinite von mises–fisher mixture modeling of whole brain
fmri data. Neural Computation, 29(10):2712–2741, 2017.
[23] H. Salehian, R. Chakraborty, E. Ofori, D. Vaillancourt,
and B. C. Vemuri. An efficient recursive estimator of the
Fréchet mean on hypersphere with applications to Medical
Image Analysis. In Math. Foundations of Computational
Anatomy, 2015.
[24] A. Schwartzman. Lognormal distributions and geometricaverages of symmetric positive definite matrices. Int. Stat.
Review, 84(3):456–486, 2016.
[25] N. Silberman, D. Hoiem, P. Kohli, and R. Fergus. Indoor
segmentation and support inference from rgbd images. In
European Conference on Computer Vision (ECCV), 2012.
[26] V. Smidl and A. Quinn. On Bayesian principal component
analysis. Computational statistics & data analysis, 51(9):4101–4123, 2007.
[27] J. Straub, T. Campbell, J. P. How, and J. W. Fisher III.
Small-variance nonparametric clustering on the hypersphere. In IEEE conference on Computer Vision and
Pattern Recognition (CVPR), 2015.
[28] 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 andStatistics, 2015.
[29] P. Thomas Fletcher. Geodesic regression and the theoryof least squares on Riemannian manifolds. International
Journal of Computer Vision (IJCV), 105(2):171–185, 11
2013.
[30] J. Traa and P. Smaragdis. Multiple speaker tracking withthe factorial von Mises-Fisher filter. In IEEE International
Workshop on Machine Learning for Signal Processing(MLSP), 2014.
[31] M. Zhang and P. Fletcher. Probabilistic Principal Geodesic
Analysis. In Advances in Neural Information Processing
Systems (NIPS) 26, pages 1178–1186, 2013.