Small-Variance Nonparametric Clustering on the Hypersphere
Julian Straub Trevor Campbell Jonathan P. How John W. Fisher III
CSAIL and LIDS, Massachusetts Institute of Technology

arXiv:1607.06407v1 [] 21 Jul 2016

{jstraub, fisher}@csail.mit.edu, {tdjc, jhow}@.mit.edu

Abstract

t0

Structural regularities in man-made environments reflect
in the distribution of their surface normals. Describing
these surface normal distributions is important in many
computer vision applications, such as scene understanding,
plane segmentation, and regularization of 3D reconstructions. Based on the small-variance limit of Bayesian nonparametric von-Mises-Fisher (vMF) mixture distributions,
we propose two new flexible and efficient k-means-like clustering algorithms for directional data such as surface normals. The first, DP-vMF-means, is a batch clustering algorithm derived from the Dirichlet process (DP) vMF mixture. Recognizing the sequential nature of data collection in
many applications, we extend this algorithm to DDP-vMFmeans, which infers temporally evolving cluster structure
from streaming data. Both algorithms naturally respect the
geometry of directional data, which lies on the unit sphere.
We demonstrate their performance on synthetic directional
data and real 3D surface normals from RGB-D sensors.
While our experiments focus on 3D data, both algorithms
generalize to high dimensional directional data such as protein backbone configurations and semantic word vectors.

t1

movement

t2

death

t3

birth
revival

Figure 1: Evolution of the data distribution on the sphere
which we model using a dependent Dirichlet process vonMises-Fisher mixture model (DDP-vMF-MM).

face normal statistics of man-made structures often necessitates the temporal integration of observations from a vast
data stream of varying cluster mixtures. Additionally, such
applications pose hard constraints on the amount of computational power available, as well as tight timing constraints.
We address these challenges by focusing on flexible
Bayesian nonparametric (BNP) Dirichlet process mixture
models (DP-MM) which describe the distribution of surface
normals in their natural space, the unit sphere in 3D. Taking
the small-variance asymptotic limit of this DP-MM of vonMises-Fisher (vMF) distributions, we obtain a fast k-meanslike algorithm, which we call DP-vMF-means, to perform
nonparametric clustering of data on the unit hypersphere.
Furthermore, we propose a novel dependent DP mixture
of vMF distributions to achieve integration of directional
data into a temporally consistent streaming model (shown
in Fig. 1). Small-variance asymptotic analysis yields the
k-means-like DDP-vMF-means algorithm. Finally, we propose a method, inspired by optimistic concurrency control,
for parallelizing the inherently sequential labeling process
of BNP-derived algorithms. This allows real-time processing of batches of 300k data-points at 30 Hz.

1. Introduction
Man-made environments and objects exhibit clear structural regularities such as planar or rounded surfaces. These
properties are evident on all scales from small objects such
as books, to medium-sized scenes like tables, rooms and
buildings and even to the organization of whole cities. Such
regularities can be captured in the statistics of surface normals that describe the local differential structure of a shape.
These statistics contain valuable information that can be
used for scene understanding, plane segmentation, or to regularize a 3D reconstruction.
Inference algorithms in fields such as robotics or augmented reality, which would benefit from the use of surface
normal statistics, are not generally provided a single batch
of data a priori. Instead, they are often provided a stream of
data batches from depth cameras. Thus, capturing the sur-

Beyond the aforementioned vision applications, directional data is ubiquitous in many other fields, including protein backbone configurations in computational biology, semantic word vectors in natural language processing, and rotations expressed as quaternions in robotics. Further, many
of those sources are observed as a stream. The proposed
algorithms directly generalize to these other data sources as
well as higher-dimensional spherical data.
1

2. Related Work
Directional distributions: A variety of distributions [30]
has been proposed in the field of directional statistics to
model data on the unit sphere. Examples are the antipodal
symmetric Bingham distribution [5], the anisotropic Kent
distribution [25], and the isotropic von-Mises-Fisher (vMF)
distribution [14]. Because of its comparative simplicity, the
vMF distribution is most commonly used.
vMF mixture models: vMF mixture models (vMF-MM)
are especially popular for modeling and inference purposes.
Banerjee et al. [3] perform Expectation Maximization (EM)
for a finite vMF mixture model to cluster text and genomic
data. This method is related to the spherical k-means (spkm)
algorithm [11], which can be obtained from a finite vMFMM by taking the infinite limit of the concentration parameter [3]. Zhong [43] extends the spherical k-means algorithm
to an online clustering framework by performing stochastic
gradient descent on the spkm objective. This approach requires the number of clusters to be known, does not allow
the creation or deletion of clusters, and heuristically weights
the contribution of old data. Gopal et al. [16] derive variational as well as collapsed Gibbs sampling inference for finite vMF-MMs, for a finite hierarchical vMF-MM and for a
finite temporally evolving vMF-MM. The vMF distribution
has also been used in BNP and hierarchical MMs. Bangert
et al. [4] formulate an infinite vMF-MM using a Dirichlet process (DP) prior, but their sampling-based inference is
known to have convergence issues and to be inefficient [23].
Reisinger et al. [34] formulate a finite, latent Dirichlet allocation (LDA) model [7] for directional data using the vMF
distribution. To the best of our knowledge there are no
other related k-means-like algorithms for directional batch
data besides the spkm algorithm. In many applications the
spherical manifold of the data is ignored and clustering is
performed in the ambient Euclidean space [19].
Surface normal modeling: In the context of object and
3D shape representation, extended Gaussian images (EGI)
have been studied [20, 21]. The EGI of a surface is the distribution of surface normals where each surface normal is
weighted by the area of the surface it represents. For computations the EGI was approximated via a tessellation of the
sphere. A crucial property of EGI is that the representation
is invariant to translation. Exploiting this property, Makadia et al. [29] extract maxima in the EGI and use the spherical FFT to compute an initial rotation estimate for pointcloud registration. Note that these maxima correspond to
the cluster centers which the proposed (D)DP-vMF-means
algorithm extracts. Furukawa et al. [15] use a tessellation of the unit sphere to extract the dominant directions
in a scene as part of a depth regularization algorithm based
on the Manhattan World (MW) assumption [10]. Tribel et
al. [42] employ EM to perform plane segmentation using

surface normals in combination with 3D locations. Assuming a finite vMF-MM and using a hierarchical clustering
approach, Hasnat et al. [17] model surface normal distributions. Straub et al. [39] propose a probabilistic model to
describe mixtures of Manhattan frames (MMF) as orthogonally coupled clusters of surface normals on the sphere. In
a similar approach, but without MMF constraints, Straub et
al. [38] introduce efficient inference for a DP-MM of clusters in distinct tangent spaces to the sphere.
Small-variance asymptotics: A principled way to reduce
the computational cost of inference, while still drawing on
the capabilities of BNP models, is to apply small-variance
analysis to the model. This technique was first applied to the
DP Gaussian mixture [26], resulting in a k-means-like algorithm with a linear regularization on the number of clusters. Since then, the technique has been extended to develop asymptotic algorithms for mixtures with general exponential family likelihoods [24], HMMs with an unknown
number of states [37], and dynamic, time-dependent mixtures [8]. However, small-variance analysis has, to date,
been limited to discrete and Euclidean spaces.

3. Von-Mises-Fisher Mixture Models
The von-Mises-Fisher (vMF) [14] distribution with
mean direction µ and concentration τ is an isotropic distribution over the D-dimensional unit hypersphere, SD−1 =
{x ∈ RD : xT x = 1}, (Fig. 2). Its density is defined as [4]
vMF(x; µ, τ ) = Z(τ ) exp(τ µT x)
−D/2

Z(τ ) = τ D/2−1 (2π)

ID/2−1 (τ )

−1

,

(1)

where µ, x ∈ SD−1 , τ > 0 and Iν is the modified Bessel
function [1] of the first kind of order ν. The conjugate prior
for µ given a fixed τ is a vMF distribution vMF(µ; µ0 , τ0 )
and the corresponding posterior given data x = {xi }N
i=1 is
p(µ|x; τ, µ0 , τ0 ) ∝ p(µ; µ0 , τ0 )

N
Y

p(xi |µ; τ )

i=1

∴ p(µ|x; τ, µ0 , τ0 ) = vMF(µ;

(2)

µ0
0
kµ0 k2 , kµ k2 ) ,

PN
where µ0 = τ0 µ0 + τ i=1 xi .
A finite mixture of K vMF distributions with known
concentration τ may be obtained by placing a Dirichlet distribution prior Dir(α) on the mixture weights π, and a vMF
prior on the mean directions µk . Data points x are assigned
to clusters via latent indicator variables z = {zi }N
i=1 .
π ∼ Dir(α), µk ∼ vMF(µ0 , τ0 ) ∀k ∈ {1, . . . , K}
zi ∼ Cat(π), xi ∼ vMF(µzi , τ ) ∀i ∈ {1, . . . , N } .

(3)

Parallel to the connection between k-means and the
Gaussian mixture model in the small-variance asymptotic

limit [26], taking τ → ∞ yields deterministic updates as
also previously noted in [3]. For completeness of the presentation we give a detailed derivation in the supplement.

4. Dirichlet Process vMF-MM
The Dirichlet process (DP) [13, 41] has been widely
used as a prior for mixture models with a countably-infinite
set of clusters [2, 4, 9, 32]. Assuming a base distribution vMF(µ; µ0 , τ0 ), the DP is an appropriate prior for a
vMF mixture with an unknown number of components and
known vMF concentration τ . Gibbs sampling inference
only differs from the finite Dirichlet vMF-MM in the label sampling step; the mixing weights π are integrated out,
resulting in the Chinese Restaurant Process (CRP) [6, 32]
p(zi = k|z−i , µ, x; τ, µ0 , τ0 )

|Ik | vMF(xi |µk ; τ ) k ≤ K
∝
α p(xi ; τ, µ0 , τ0 )
k = K + 1,

(4)

where Ik is the set of data indices assigned to cluster k.
Note that the DP concentration parameter, α > 0, influences the likelihood of adding a new clusters. The conjugate prior for µk yields p(xi ; τ, µ0 , τ0 ) via marginalization:
Z
p(xi ; τ, µ0 , τ0 ) = vMF(xi |µk ; τ ) vMF(µk ; µ0 , τ0 ) dµk
=

Z(τ )Z(τ0 )
.
Z(||τ xi + τ0 µ0 ||2 )

(5)

4.1. DP-vMF-means
In this section, we provide a small-variance asymptotic
analysis of the label and parameter update steps of the Gibbs
sampling algorithm for the DP vMF mixture, yielding deterministic updates. As with k-means, the label assignments are computed sequentially for all datapoints before
the means are updated, and the process is iterated until convergence. Pseudocode can be found in the supplement.
Label Update: To derive a hyperspherical analog to DPmeans [26], consider the limit of the label sampling step (4)
as τ → ∞. The normalizer Z(||τ xi + τ0 µ0 ||2 ) approaches
τ (∗)

Z(||τ xi + τ0 µ0 ||2 ) → exp(−τ ) ,

(6)

where overscript τ (∗) denotes proportionality up to a finite
power of τ , and where we have used the fact that as τ → ∞,
the modified Bessel function of the first kind satisfies [1]

 
1
exp(τ )
τ (∗)
ID/2−1 (τ ) = √
1−O
→ exp (τ ) . (7)
τ
2πτ
To achieve a nontrivial result, the asymptotic behavior of
Z(τ ) must be matched by α, so let α = exp(λτ ) to obtain
αZ(τ )Z(τ0 )
τ (∗)
= Z(τ ) exp(τ (λ + 1)) .
Z(||τ xi + τ0 µ0 ||2 )

(8)

Therefore, as τ → ∞, the label sampling step becomes
lim p(zi = k|z−i , µ, x; τ, µ0 , τ0 )

τ (xT µk −λ−1)

 PK |Ik |e τ (xiT µ −λ−1)
i j
+c(τ )
j=1 |Ij |e
= lim
c(τ )
τ →∞
 PK
T
τ (x µj −λ−1)

τ →∞

j=1

|Ij |e

i

+c(τ )

k≤K

(9)

k = K + 1,

where we have used that the normalizers Z(τ ) of
τ (∗)

vMF(xi |µk ; τ ) and Eq. (8) cancel and c(τ ) = 1. Thus, as
τ → ∞, sampling from p(zi |z−i , µ, x; τ, µ0 , τ0 ) is equivalent to the following assignment rule:
 T
xi µk k ≤ K
zi = arg max
(10)
λ + 1 k = K + 1.
k∈{1,...,K+1}
Since −1 ≤ xTi µk ≤ 1, the parameter λ can be restricted
to the set λ ∈ [−2, 0] without loss of generality. Intuitively
λ defines the maximum angular spread φλ of clusters about
their mean direction, via λ = cos(φλ ) − 1. Note, that upon
assigning a datapoint to a new cluster, i.e. zi = K + 1, the
mean of that cluster is initialized to µK+1 = xi . Finally, if
an observation xi is the last one in its cluster, the cluster is
removed prior to finding the new label for xi using Eq. (10).
Parameter Update: Taking τ → ∞ in the parameter posterior for cluster k from Eq. (2) causes τ0 and µ0 to become
negligible. Hence the parameter update becomes:
P
xi
µk = P i∈Ik
∀k ∈ {1, . . . , K} .
(11)
k i∈Ik xi k2
Objective Function: From Eq. (10) we can see that assigning a datapoint xi to cluster k provides a score of
xTi µk , whereas adding a new cluster provides a score of
λ+1−xTi µK+1 = λ, since new mean directions are initialized directly to µK+1 = xi . Hence, the objective function
that DP-vMF-means maximizes is
JDP-vMF =

K X
X

xTi µk + λK .

(12)

k=1 i∈Ik

5. Dependent Dirichlet Process vMF-MM
Suppose now that, in addition to an unknown number K
of components, the vMF mixture undergoes temporal evolution in discrete timesteps t ∈ N (Fig. 1): mixture components can move, be destroyed, and new ones can be created
at each timestep. For such a scenario, the dependent Dirichlet process (DDP) [28, 27, 8] is an appropriate prior over the
mixture components and weights. Using intermediate auxiliary DPs F0 and F1 , the DDP constructs a Markov chain
of DPs Gt , where Gt+1 is sampled from Gt as follows:
1. (Death) For each atom θ in Gt , sample from
Bernoulli(q). If the result is 1, add θ to F0 .

2. (Motion) Replace each θ in F0 with θ0 ∼ T (θ0 |θ).
3. (Birth) Sample a DP F1 ∼ DP(α, H). Let Gt+1 be a
random convex combination of F0 and F1 .

walk µk0 → µk1 → · · · → µk∆tk = µk with initial distribution µk0 ∼ vMF(µk0 ; mk , τk ). Therefore, the intermek
diate mean directions {µkn }∆t
n=1 must be marginalized out
when computing p(xi ; mk , τk ):

There are four parameters in this model: α > 0 and H(·),
the concentration parameter and base measure of the innovation process; q ∈ (0, 1), the Bernoulli cluster survival
probability; and finally T (·|·), the random walk transition
distribution. In the present work, both the base and random transition distributions are von-Mises-Fisher: H(µ) =
vMF (µ; µ0 , τ0 ), and T (µ|ν) = vMF (µ; ν, ξ).
Suppose at timestep t, a new batch of data x is observed.
Then Gibbs sampling posterior inference for the DDP mixture, as in the previous sections, iteratively samples labels and parameters. Let the set of tracked mean directions from previous timesteps be {µk0 }K
k=1 , where µk0 ∼
vMF(mk , τk ) and ∆tk denotes the number of timesteps
since cluster k was last instantiated (i.e. when it last had
data assigned to it).1 Then the label sampling distribution is

The integration in (15) cannot be computed in closed
form; however, the value of the integral is only of interest
in the limit as τ → ∞. Therefore, Theorem 1, an extension
of Laplace’s approximation to general differentiable manifolds, may be used to obtain an exact formula.

p(zi = k|µ, z−i , x)
 1−qt
k =K +1
 α 1−q p(xi ; µ0 , τ0 )
∝ (ck + |Ik |) vMF(xi |µk ; τ ) k ≤ K, |Ik | > 0
 ∆tk
q ck p(xi ; mk , τk )
k ≤ K, |Ik | = 0 ,

Theorem 1 (Manifold Laplace Approximation). Suppose
M ⊂ Rn is a bounded m-dimensional differentiable manifold and f : Rn → R is a smooth function on M . Further, suppose f has a unique global maximum on M , x? =
arg maxx∈M f (x). Then

(13)

where ck is the number of observations assigned to cluster
k in past timesteps. The parameter sampling distribution is

p(xi ; mk , τk )
R
R p(xi |µk∆tk ; τ ) · p(µ0 ; mk , τk )
Q∆tk
=
···
· n=1
p(µkn |µk(n−1) ; ξ)
µk0 ,...,µk∆tk
R
R
∆tk
= Z(τ )Z(βτ ) Z(wk τ )
···
exp (τ f )
f = xTi µk∆tk + β

∆tk
X

µTkn µk(n−1) + wk µTk0 mk .

n=1

R
lim

τ →∞

M


2π m
τ

(
µ0
vMF(µk ; kµ0kk2 , kµ0k k2 ) ck = 0
k
p(µk |µ−k , z, x) ∝
p(µk |x, z; mk , τk )
ck > 0 ,

(15)

µk0 ,...,µk∆tk

eτ f (x)
1

|det U T ∇2 f (x? )U | 2 eτ f (x? )

= 1,

(16)

(14)

where U ∈ Rn×m is a matrix whose columns are an orthonormal basis for the tangent space of M at x? .

P
where µ0k = τ0 µ0 + τ i∈Ik xi , and p(µk |x, z; mk , τk ) is
the distribution over the current cluster k mean direction µk
given the assigned data and the old mean direction µk0 .

Proof. See the supplementary material. The general technique of this proof is to transform coordinates between
the manifold and its tangent plane using the exponential
map [12], and then apply the multidimensional Laplace approximation in the transformed Euclidean space.

5.1. DDP-vMF-means
In the following, we analyze the small-variance asymptotics of the DDP-vMF mixture model. We first derive the
label assignment rules, followed by the parameter updates.
Label Update: First, let α = exp(λτ ), q = exp(Qτ ),
ξ = exp(βτ ), and τk = τ wk , with λ ∈ [−2, 0] as before,
t
Q ≤ 0, and β, wk ≥ 0. Note that limτ →∞ 1−q
1−q = 1, and
thus the asymptotics of the label assignment probability for
current and new clusters is the same as in Section 4.1.
Hence, we focus on the assignment of a datapoint to a
previously observed, but currently not instantiated, cluster
k. During the ∆tk timesteps since cluster k was last observed, the mean direction µk underwent a random vMF
1 Note

that all quantities (x, µk0 , mk , τk , ck , nk , etc.) are now timevarying. This dependence is not shown in the notation for brevity, and all
quantities are assumed to be shown for the current timestep t.

N
Corollary 1. Given a smooth function f : RD → R, with
a unique global maximum over the N -product of (D − 1)N
spheres x? ∈ SD−1 ,
Z
?
τ (∗)
eτ f (x) → eτ f (x ) .
(17)
(SD−1 )N

Proof. This is Theorem 1 applied to SD−1

N

.

Using Corollary 1 and the limiting approximation of the
modified Bessel function (7) in equation (15) yields the following asymptotic behavior as τ → ∞:
τ (∗)

p(xi ; mk , τk ) → exp (τ (f ? − 1 − β∆tk − wk )) .
(18)

ζ

θ

η

(a) τ = 100

mk
µk0
φ
··
φ

·
·
·
µk∆tk
x̄k
||x̄k ||2

(b) τ = 1

Figure 2: Left: 2D vMF distributions. Right: geometry of
the maximum likelihood setting of µk0 , µk1 , . . . , µk∆tk for
the transition distribution.
The only remaining unknown in the asymptotic expression, f ? , can be found via constrained optimization
max

∆tk
{µkn }n=1

xTi µk∆tk + β

∆tk
X

µTkn µk(n−1) + wk µTk0 mk

n=1

(19)

s.t. µTkn µkn = 1 ∀n ∈ {0, . . . , ∆tk } .
The optimization (19) has a closed-form solution:
wk mk + βµk1
||wk mk + βµk1 ||2
µk(n+1) + µk(n−1)
=
||µk(n+1) + µk(n−1) ||2
xi + βµk(∆tk −1)
.
=
||xi + βµk(∆tk −1) ||2

µk0 =
µkn
µk∆tk

∀n∈{1,...,∆tk −1}

(20)

rule zi = arg maxk Jk , where


λT+ 1


µk xi


∆tk β(cos(φ? ) − 1)
Jk =


 +wk (cos(θ? ) − 1) 



+ cos(η ? ) + ∆tk Q

Parameter Update: The parameter update rule for DDPvMF-means comes from the asymptotic behavior of (14) as
τ → ∞. The analysis for any new cluster is the same as
that in Section 4.1, so our focus is again on the transitioned
mean direction posterior p(µk∆tk |x, z; mk , τk ) (recall that
µk = µk∆tk in the definition of the random vMF walk).
This distribution can be expanded, similarly to (15), as:
p(µk∆tk |x, z; mk , τk )
R
R
p(x|µk∆tk ; τ ) · p(µ0 ; mk , τk )
Q∆tk
=
···
µk0 ,...,µk(∆tk −1) ·
n=1 p(µkn |µk(n−1) ; ξ)
R
R
|Ik |
∆tk
= Z(τ ) Z(βτ ) Z(wk τ )
···
exp (τ f )

?

?

ζ = θ + ∆tk φ + η =

arccos(mTk xi ) ,

X

xTi µk∆tk + β

∆tk
X

µTkn µk(n−1) + wk µTk0 mk .

n=1

i∈Ik

P
Define x̄k = i∈Ik xi . Once again, applying Corollary 1,
the limit τ → ∞ removes the integrals over the marginalized mean directions. However, in contrast to the label assignment update, µk∆tk is not marginalized out. Therefore,
an additional maximization with respect to µk∆tk to find the
concentration point of the posterior yields
τ (∗)

(21)

p(µk∆tk |x, z; mk , τk ) →

(26)

exp(τ (f ? − |Ik | − ∆tk β − wk ))

where ζ is the full angle between xi and mk . Since (21)
cannot be solved in closed-form, Newton’s method is used
to compute φ? , θ? , and η ? , which in turn determines f ? :
f ? = wk cos(θ? ) + β∆tk cos(φ? ) + cos(η ? ) .

(25)

µk0 ,...,µk(∆tk −1)

These ternary relationships enforce that the optimal vMF
mean directions along the random walk lie on the geodesic
between mk and xi . Therefore, this walk can be described
geometrically by three angles, as shown in Fig. 2 (with
x̄k = xi ): the angle φ between consecutive µkn , the angle
η between xi and µk∆tk , and the angle θ between mk and
µk0 . Given these definitions, standard trigonometry yields
a set of three equations in φ, η, and θ:

?

(24)
k ≤ K, |Ik | = 0 .

Note that if Q∆tk < λ, cluster k can be removed permanently as it will never be revived again. Furthermore, for
any Q ≤ λ all clusters are removed after each timestep, and
the algorithm reduces to DP-vMF-means.

f=

wk sin(θ? ) = β sin(φ? ) = sin(η ? )

k =K +1
k ≤ K, |Ik | > 0

(22)

Returning to (18), the transition asymptotics are


τ wk (cos(θ? ) − 1)
τ (∗)
p(xi ; mk ; τk ) → exp +τ β∆tk (cos(φ? ) − 1) . (23)
+τ (cos(η ? ) − 1)
Substituting this into Eq. (13) with the earlier definition q =
exp(τ Q), and taking the limit τ → ∞ yields the assignment

?

?

?

?

f = wk cos(θ ) + β∆tk cos(φ ) + ||x̄k ||2 cos(η ) .
Analyzing the geometry of the geodesic between x̄k /kx̄k k2
and mk (Fig. 2) there exist φ? , θ? and η ? such that
wk sin(θ? ) = β sin(φ? ) = ||x̄k ||2 sin(η ? )
ζ = θ? + ∆tk φ? + η ? = arccos(mTk kx̄x̄kkk2 ) ,

(27)

which can be solved via Newton’s method. Given the solution, µk can be obtained by rotating kx̄x̄kkk2 by angle η ? on
the geodesic shown in Fig. 2 towards mk ,
µk = R(η ? ) ||x̄x̄kk||2 .

(28)

φλ [deg]
10 15 20 25 30 35 40 45 1.00
KT =30

φλ [deg]
10 15 20 25 30 35 40 45

0.95
silhouette

DP-vMF-means
spkm

NMI

0.90

K

90
80
70
60
50
40
30
20
10
0

0.85

0.80

60

50

40

30

20

10 0.75

60

50

40

30

20

10

1.0
0.8
0.6
0.4
0.2
0.0
0.2
0.4
0.6

φλ [deg]
10 15 20 25 30 35 40 45

60

50

40

30

20

10

80

90

φλ [deg]
100 110 120 130 140

7

DP-vMF-means
spkm
6
5
4

K
K
K
Figure 3: Comparison of the spkm and the DP-vMF-means clustering algorithms on synthetic spherical data with KT = 30
clusters. Note DP-vMF-means’ higher maximum normalized mutual information (NMI) as well as silhouette score.

6. Optimistic Iterated Restarts (OIR)
In our implementation of the algorithm we pay special
attention to speed and parallel execution to enable real-time
performance for streaming RGB-D data.
Observe that the main bottleneck of DP-based hard clustering algorithms, such as the proposed (D)DP-vMF-means,
DP-means [26] or Dynamic means [8], is the inherently sequential assignment of labels: due to the creation of new
clusters, the label assignments depend on all previous assignments. While this is a key feature of the streaming clustering algorithms, it poses a computational hindrance. We
address this issue with an optimistic parallel label assignment procedure inspired by techniques for database concurrency control [33].
First, we compute assignments in parallel (e.g. on a
GPU). If all datapoints were assigned only to instantiated
clusters, we output the labeling. Otherwise, we find the
lowest observation id i that modified the number of clusters,
apply the modification, and recompute the assignments for
all observations i0 > i in parallel. Thus, per data-batch,
DP-vMF-means restarts once per new cluster, while DDPvMF-means restarts once for each new or revived cluster.

7. Results
7.1. Evaluation of the DP-vMF-means Algorithm
Synthetic Data: First, we evaluate the
behavior of the DP-vMF-means algorithm
in comparison to its parametric cousin, the
spkm algorithm, on synthetic 3D spherical data sampled from KT = 30 true
vMF distributions. All evaluation results
Synthetic Data
are shown as the mean and standard de-

900
800
700
600
500
400
300
200
100
02 3 4 5 6 7 8 9 1011
number of clusters

0.85
0.80
0.75
0.70
0.65
0.60
0.55
0.50

silhouette

Weight Update: After the iteration of label and parameter updates has converged, the weight wk must be updated for all clusters to reflect the new uncertainty in the
mean direction of cluster k. This can be done by examining
(26): Since at the maximum of a vMF(µ; mk , wk τ ) density,
exp(τ wk mTk µ) = exp(τ wk ), wk is updated to f ? .

8

K

3

2

Figure 4: Histogram over the number of clusters found by
DP-vMF-means (left) and silhouette values for DP-vMFmeans and spkm (right) across the whole NYU dataset [31].
viation over 50 runs. The left plot of Fig. 3 depicts the
inferred number of clusters K on the horizontal axis as a
function of the respective parameters of the two algorithms:
the number of clusters K for spkm and the parameter φλ
for DP-vMF-means (recall that φλ = cos−1 (λ + 1) as defined in Sec. 4.1). This figure demonstrates the ability of the
DP-vMF-means algorithm to discover the correct number of
clusters KT , and the relative insensitivity of the discovered
number of clusters with respect to its parameter φλ .
The middle and right hand plots show two measures
for clustering quality. The normalized Mutual Information (NMI) [40], depicted in the middle, is computed using the true labels. DP-vMF-means achieves an almost perfect NMI of 0.99, while spkm only reaches 0.94 NMI even
with K = KT . The slightly superior performance of DPvMF-means stems from its enhanced ability to avoid local
optima due to the way labels are initialized: while spkm is
forced to initialize K cluster parameters, DP-vMF-means
starts with an empty set and adds clusters on the fly as more
data are labelled. The NMI results are corroborated by the
silhouette score [36], shown to the right in Fig. 3. The silhouette score is an internal measure for clustering quality
that can be computed without knowledge of the true clustering, and is used to tune parametric clustering algorithms.
With a maximum of 0.92 DP-vMF-means reaches a close to
perfect silhouette score, indicating well-separated, concentrated clusters. Again, spkm does not reach the same clus-

RGB
DP-vMF-means
K=4
spkm K = 5

Figure 5: Directional segmentation of scenes from the NYU v2 RGB-D dataset [31] as implied by surface normal clusters.
The complexity of the scenes increases from left to right as can be seen from the RGB images. The second row shows the
clustering inferred by DP-vMF-Means while the third and fourth show the spherical k-means results. Black denotes missing
data due to sensor limitations. Note that DP-vMF-means adapts the number of clusters to the complexity of the scene.
tering performance even for K = KT for the same aforementioned reasons.
For additional reference we ran two sampling-based inference algorithms for the DP-vMF-MM model. The CRPbased inference of [4] was aborted after running for two
days without convergence. This inefficiency has been noted
previously [23]. A more efficient alternative is the finite
Dirichlet process (FSD) approximation [22] to the DPvMF-MM. The inference can be parallelized and yielded
results within minutes. After 5000 iterations the sampler
converged to an incorrect number of 23.0 ± 1.67 clusters
on average over 10 sampler runs. These convergence issues
have also been noted before [9]. The last sample of the different runs achieves an average NMI of 0.94 ± 0.01 and a
silhouette score of 0.74 ± 0.04.
NYU v2 depth dataset: In this experiment, the DPvMF-means and spkm algorithms were compared on the
NYU v2 RGB-D dataset [31]. Surface normals were extracted from the depth images [19] and preprocessed with
total variation smoothing [35]. We quantify the clustering quality in terms of the average silhouette score over the
clusterings of the 1449 scenes of the NYU v2 depth dataset.
Since we do not possess the true scene labeling, we use the
silhouette quality metric as a proxy for the NMI metric; this
was motivated by results of the synthetic experiment.
Across the whole NYU v2 dataset, the DP-vMF-means
algorithm achieves the highest average silhouette score of
0.75 for φ?λ = 100◦ as depicted in Fig. 4. The histogram
over the number of inferred clusters by DP-vMF-means for
φ?λ indicates the varying complexity of the scenes ranging

from three to eleven. The clear peak at K = 4 coincides
with the highest silhouette score for spkm (0.73) and explains the only slightly lower silhouette score of spkm: most
scenes in the dataset exhibit four primary directions.
Figure 5 shows a qualitative comparison of the scene
segmentation implied by the clustering of surface normals.
In comparison to spkm, the DP-vMF-means clustering results show the ability of the algorithm to adapt the number of clusters to the scene at hand. If the right number of
clusters is selected for the spkm clustering, the results have
similar quality; however, the number of clusters is generally
not known a priori and varies across scenes. This demonstrates two major advantages of DP-vMF-means over spkm:
(1) DP-vMF-means is less sensitive to the parameter setting
(see Fig, 3, left) and (2) it is easier to choose φλ than K
since it intuitively corresponds to the maximum angular radius of a cluster, which can be gauged from the type of data
and its noise characteristics. For this experiment φλ = 100◦
is justified by the typical Manhattan structure [10] of the indoor environment plus 10◦ to account for sensor noise.

7.2. Evaluation of the DDP-vMF-means Algorithm
Real-time Directional Segmentation: In fields such as
mobile robotics or augmented reality, it is uncommon to observe just a single RGB-D frame of a scene; more typically,
the sensor will observe a temporal sequence of frames. The
following experiment demonstrates the temporally consistent clustering capability of the DDP-vMF-means algorithm
on surface normals extracted from a sequence of depth images recorded in an indoor environment. Each frame is pre-

segmentation

DDP-vMF-means

segmentation

DP-vMF-means

segmentation

sDP-vMF-means

cluster shares

60
10
60
10

30

120

210

frame number 330

420

510

600

Figure 6: Clustering of a surface normal stream recorded when walking a 90◦ turn in an office environment. We depict
key-frames color-coded with the implied surface-normal clustering for three clustering algorithms. The plots in the second
and third row depict the percentage of normals associated to the respective cluster for DDP-vMF-means and sequential DPvMF-means. Note that only the clustering obtained via the DDP-vMF-means algorithm is consistent across the whole run.

processed in 11ms using edge-preserving smoothing with a
hybrid CPU-GPU guided filter [18].
We compare against the ad-hoc approaches of clustering on a frame-by-frame basis using DP-vMF-means, both
with and without initializing the algorithm from the previous frame’s clusters. The former is referred to as sequential DP-vMF-means (sDP-vMF-means). sDP-vMF-means
achieves a greedy frame-to-frame label consistency, but, unlike DDP-vMF-means, it cannot reinstantiate previous clusters after multiframe lapses. Motivated by the DP-vMFmeans evaluation, all algorithms were run with φλ = 100◦ .
λ
For DDP-vMF-means β = 105 and Q = 400
.
The differences in labeling consistency can be observed
in rows two and three of Fig. 6, which shows the percentage
of normals associated with a specific cluster. While DDPvMF-means is temporally consistent and reinstantiates the
lime-green and red clusters, observed in the first half of the
run, DP-vMF-means erroneously creates new clusters. We
do not depict the percentages of surface normals associated
with the clusters for the batch DP-vMF-means algorithm,
since the there is no label consistency between time-steps
as can be observed in the last row of Fig. 6.
The average run-time per frame was 28.4 ms for batch
DP-vMF-means, 12.8 ms for sDP-vMF-means, 20.4 ms for
DDP-vMF-means, and 13.6 ms for spkm with K = 5. The
increased running time of batch DP-vMF-means is a result
of clustering each batch of surface normals in isolation; OIR
label assignment needs several restarts to assign labels to all
surface normals. By initializing the clusters from a previ-

ous frame, sDP-vMF-means only incurs labeling restarts if
a new cluster is observed, and hence has significantly lower
run time. DDP-vMF-means is slightly slower than sDPvMF-means since it is keeping track of both observed and
unobserved clusters.

8. Conclusion
Taking the small-variance asymptotic limit of the
Bayesian nonparametric DP-vMF and DDP-vMF mixture
models, we have derived two novel spherical kmeans-like
algorithms for efficient batch and streaming clustering on
the unit hypersphere. The performance and flexibility of
DP-vMF-means was demonstrated on both synthetic data
and the NYU v2 RGB-D dataset. For DDP-vMF-means,
Optimistic Iterated Restarts (OIR) parallelized label assignments, enable real-time temorally consistent clustering of
batches of 300k surface normals collected at 30 Hz from a
RGB-D camera.
We envision a large number of potential applications for
the presented algorithms in computer vision and in other
realms where directional data is encountered. Implementations are available at http://people.csail.mit.
edu/jstraub/.
Acknowledgements
This work was partially supported by ONR MURI N0001411-1-0688 and ARO MURI W911NF-11-1-0391.

References
[1] M. Abramowitz and I. Stegun, editors. Handbook of Mathematical Functions. Dover Books on Mathematics. Dover
Publications, 1965.
[2] C. Antoniak. Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. The Annals of
Statistics, 1974.
[3] A. Banerjee, I. S. Dhillon, J. Ghosh, S. Sra, and G. Ridgeway.
Clustering on the unit hypersphere using von Mises-Fisher
distributions. JMLR, 6(9), 2005.
[4] M. Bangert, P. Hennig, and U. Oelfke. Using an infinite von
Mises-Fisher mixture model to cluster treatment beam directions in external radiation therapy. In ICMLA, 2010.
[5] C. Bingham. An antipodally symmetric distribution on the
sphere. The Annals of Statistics, 2(6):1201–1225, 1974.
[6] D. Blackwell and J. B. MacQueen. Ferguson distributions
via pólya urn schemes. The Annals of Statistics, 1973.
[7] D. M. Blei, A. Y. Ng, and M. I. Jordan. Latent Dirichlet
allocation. JMLR, 3:993–1022, 2003.
[8] T. Campbell, M. Liu, B. Kulis, J. P. How, and L. Carin. Dynamic clustering via asymptotics of the dependent Dirichlet
process mixture. In NIPS, 2013.
[9] J. Chang and J. W. Fisher III. Parallel sampling of dp mixture
models using sub-clusters splits. In NIPS, 2013.
[10] J. M. Coughlan and A. L. Yuille. Manhattan world: Compass direction from a single image by Bayesian inference. In
ICCV, 1999.
[11] I. S. Dhillon and D. S. Modha. Concept decompositions for
large sparse text data using clustering. Machine learning,
42(1-2):143–175, 2001.
[12] M. P. do Carmo. Riemannian Geometry. Birkhäuser Verlag,
Boston, MA, 1992.
[13] T. Ferguson. A Bayesian analysis of some nonparametric
problems. The Annals of Statistics, 1973.
[14] N. I. Fisher. Statistical Analysis of Circular Data. Cambridge
University Press, 1995.
[15] Y. Furukawa, B. Curless, S. M. Seitz, and R. Szeliski.
Manhattan-world stereo. In CVPR, 2009.
[16] S. Gopal and Y. Yang. von Mises-Fisher clustering models.
In ICML, 2014.
[17] M. A. Hasnat, O. Alata, and A. Trémeau. Hierarchical 3-d
von Mises-Fisher mixture model. In Workshop on Divergences and Divergence Learning, ICML, 2013.
[18] K. He, J. Sun, and X. Tang. Guided image filtering. In ECCV,
2010.
[19] D. Holz, S. Holzer, and R. B. Rusu. Real-Time Plane Segmentation using RGB-D Cameras. In Proceedings of the
RoboCup Symposium, 2011.
[20] B. K. P. Horn. Extended Gaussian images. Proceedings of
the IEEE, 72(12):1671–1686, 1984.
[21] K. Ikeuchi. Recognition of 3-d objects using the extended
Gaussian image. In In IJCAI Conference, 1981.
[22] H. Ishwaran and M. Zarepour. Exact and approximate sum
representations for the Dirichlet process. Canadian Journal
of Statistics, 30(2):269–283, 2002.

[23] S. Jain and R. Neal. A split-merge Markov chain Monte
Carlo procedure for the Dirichlet process mixture model.
Journal of Computational and Graphical Statistics, 13:158–
182, 2000.
[24] K. Jiang, B. Kulis, and M. Jordan. Small-variance asymptotics for exponential family Dirichlet process mixture models. In NIPS, 2012.
[25] J. T. Kent. The Fisher-Bingham distribution on the sphere.
Journal of the Royal Statistical Society, pages 71–80, 1982.
[26] B. Kulis and M. I. Jordan. Revisiting k-means: New algorithms via Bayesian nonparametrics. In ICML, 2012.
[27] D. Lin, E. Grimson, and J. Fisher. Construction of dependent
Dirichlet processes based on Poisson processes. NIPS, 2010.
[28] S. N. MacEachern. Dependent nonparametric processes. In
ASA Proceedings of the Section on Bayesian Statistical Science, 1999.
[29] A. Makadia, A. Patterson, and K. Daniilidis. Fully automatic
registration of 3D point clouds. In CVPR, 2006.
[30] K. V. Mardia and P. E. Jupp. Directional statistics, volume
494. John Wiley & Sons, 2009.
[31] P. K. Nathan Silberman, Derek Hoiem and R. Fergus. Indoor
segmentation and support inference from RGBD images. In
ECCV, 2012.
[32] R. Neal. Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 9(2):249–265, 2000.
[33] X. Pan, J. E. Gonzalez, S. Jegelka, T. Broderick, and M. Jordan. Optimistic concurrency control for distributed unsupervised learning. In NIPS, 2013.
[34] J. Reisinger, A. Waters, B. Silverthorn, and R. J. Mooney.
Spherical topic models. In ICML, 2010.
[35] G. Rosman, Y. Wang, X.-C. Tai, R. Kimmel, and A. M.
Bruckstein. Fast regularization of matrix-valued images. In
ECCV, 2012.
[36] P. J. Rousseeuw. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics, 20:53–65, 1987.
[37] A. Roychowdhury, K. Jiang, and B. Kulis. Small-variance
asymptotics for hidden Markov models. In NIPS, 2013.
[38] J. Straub, J. Chang, O. Freifeld, and J. W. Fisher III. A
Dirichlet process mixture model for spherical data. In AISTATS, 2015.
[39] J. Straub, G. Rosman, O. Freifeld, J. J. Leonard, and J. W.
Fisher III. A mixture of Manhattan frames: Beyond the Manhattan world. In CVPR, 2014.
[40] A. Strehl and J. Ghosh. Cluster ensembles—a knowledge
reuse framework for combining multiple partitions. JMLR,
3:583–617, 2003.
[41] Y. W. Teh. Dirichlet processes. In Encyclopedia of Machine
Learning. Springer, 2010.
[42] R. Triebel, W. Burgard, and F. Dellaert. Using hierarchical
EM to extract planes from 3D range scans. In ICRA, 2005.
[43] S. Zhong. Efficient online spherical k-means clustering. In
IJCNN, volume 5, pages 3180–3185. 2005.

