A Sub-Quadratic Exact Medoid Algorithm

arXiv:1605.06950v4 [stat.ML] 12 Apr 2017

James Newling
Idiap Research Institute & EPFL

Abstract
We present a new algorithm trimed for obtaining the medoid of a set, that is the element of the set which minimises the mean
distance to all other elements. The algorithm
is shown to have, under certain assumptions,
3
expected run time O(N 2 ) in Rd where N is
the set size, making it the first sub-quadratic
exact medoid algorithm for d > 1. Experiments show that it performs very well on spatial network data, frequently requiring two
orders of magnitude fewer distance calculations than state-of-the-art approximate algorithms. As an application, we show how
trimed can be used as a component in an
accelerated K-medoids algorithm, and then
how it can be relaxed to obtain further computational gains with only a minor loss in
cluster quality.

1

Introduction

A popular measure of the centrality of an element of
a set is its mean distance to all other elements. In
network analysis, this measure is referred to as closeness centrality, we will refer to it as energy. Given
a set S = {x(1), . . . , x(N )} the energy of element
i ∈ {1, . . . , N } is thus given by,
E(i) =

1
N

X

dist(x(i), x(j)).

j∈{1,...,N }

An element in S with minimum energy is referred to
as a 1-median or a medoid. Without loss of generality, we will assume that S contains a unique medoid.
The problem of determining the medoid of a set arises
in the contexts of clustering, operations research, and
Proceedings of the 20th International Conference on Artificial Intelligence and Statistics (AISTATS) 2017, Fort Lauderdale, Florida, USA. JMLR: W&CP volume 54. Copyright 2017 by the author(s).

François Fleuret
Idiap Research Institute & EPFL

network analysis. In clustering, the Voronoi iteration
K-medoids algorithm (Hastie et al., 2001; Park and
Jun, 2009) requires determining the medoid of each of
K clusters at each iteration. In operations research,
the facility location problem requires placing one or
several facilities so as to minimise the cost of connecting to clients. In network analysis, the medoid may
represent an influential person in a social network, or
the most central station in a rail network.
1.1

Medoid algorithms and our contribution

A simple algorithm for obtaining the medoid of a set
of N elements computes the energy of all elements and
selects the one with minimum energy, requiring Θ(N 2 )
time. In certain settings Θ(N ) algorithms exist, such
as in 1-D where the problem is solved by Quickselect (Hoare, 1961), and more generally on trees. However, no general purpose o(N 2 ) algorithm exists. An
example illustrating the impossibility of such an algorithm is presented in Supplementary Material B (SMA). Related to finding the medoid of a set is finding
the geometric median, which in vector spaces is defined as the point in the vector space with minimum
energy. The relationship between the two problems is
discussed in §2.1.
Much work has been done to develop approximate algorithms in the context of network analysis. The RAND
algorithm of Eppstein and Wang (2004) can be used to
estimate the energy of all nodes in a graph. The accuracy of RAND depends on the diameter of the network,
which motivated Cohen et al. (2014) to use pivoting
to make RAND more effective for large diameter networks. The work most closely related to ours is that
of Okamoto et al. (2008), where RAND is adapted to
the task of finding the k lowest energy nodes, k = 1
corresponding to the medoid problem. The resulting
TOPRANK algorithm of Okamoto et al. (2008) has run
time Õ(N 5/3 ) under certain assumptions, and returns
the medoid with probability 1 − O(1/N ), that is with
high probability (w.h.p.). Note that only their run time
result requires any assumption, obtaining the medoid
w.h.p. is guaranteed. TOPRANK is discussed in §2.2.
In this paper we present an algorithm which has ex-

A Sub-Quadratic Exact Medoid Algorithm

pected run time O(N 3/2 ) under certain assumptions
and always returns the medoid. In other words, we
present an exact medoid algorithm with improved
complexity over the state-of-the-art approximate algorithm, TOPRANK. We show through experiments that
the new algorithm works well for low-dimensional data
in Rd and for spatial network data. Our new medoid
algorithm, which we call trimed, uses the triangle inequality to quickly eliminate elements which cannot be
the medoid. The O(N 3/2 ) run time follows from the
surprising result that all but O(N 1/2 ) elements can be
eliminated in this way.
The complexity bound on expected run time which we
derive contains a term which grows exponentially in
dimension d, and experiments show that in very high
dimensions trimed often ends up computing O(N 2 )
distances.
1.2

K-medoids algorithms and our
contribution

It should be noted that algorithms other than KMEDS
have been proposed for finding approximate solutions
to the K-medoids problem, and have been shown to be
very effective in Newling and Fleuret (2016b). These
include PAM and CLARA of Kaufman and Rousseeuw
(1990), and CLARANS of Ng et al. (2005). In this paper we do not compare cluster qualities of previous
algorithms, but focus on accelerating the lloyd equivalent for K-medoids as a test setting for our medoid
algorithm trimed.

2
2.1

Previous works
A related problem: the geometric median

A problem closely related to the medoid problem is the
geometric median problem. In the vector space Rd the
geometric median, assuming it is unique, is defined as,


X
g(S) = arg min 
kv − yk .
(1)
v∈V

The K-medoids problem is to partition a set into K
clusters, so as to minimise the sum over elements of
dissimilarites with their nearest medoids. That is, to
choose M = {m(1), . . . , m(K)} ⊂ {1, . . . , N } to minimise,
L(M) =

N
X
i=1

min

diss(x(i), x(m(k))).

k∈{1,...,K}

We focus on the special case where the dissimilarity
is a distance (diss = dist), which is still more general
than K-means which only applies to vector spaces. Kmedoids is used in bioinformatics where elements are
genetic sequences or gene expression levels (Chipman
et al., 2003) and has been applied to clustering on
graphs (Rattigan et al., 2007). In machine vision, Kmedoids is often preferred, as a medoid is more easily
interpretable than a mean (Frahm et al., 2010).
The K-medoids problem is NP-hard, but there exist
approximation algorithms. The Voronoi iteration algorithm, appearing in Hastie et al. (2001) and later
in Park and Jun (2009), consists of alternating between updating medoids and assignments, much in the
same way as Lloyd’s algorithm works for the K-means
problem. We will refer to it as KMEDS, and to Lloyd’s
K-means algorithm as lloyd.
One significant difference between KMEDS and lloyd is
that the computation of a medoid is quadratic in the
number of elements per cluster whereas the computation of a mean is linear. By incorporating our new
medoid algorithm into KMEDS, we break the quadratic
dependency of KMEDS, bringing it closer in performance
to lloyd. We also show how ideas for accelerating
lloyd presented in Elkan (2003) can be used in KMEDS.

y∈S

While the medoid of a set is defined in any space with a
distance measure, the geometric median is specific to
vector spaces, where addition and scalar multiplication are defined. The convexity of the objective function being minimised in (1) has enabled the development of fast algorithms. In particular, Cohen et al.
(2016) present an algorithm which obtains an estimate
for the geometric median with relative error 1 + O()
with complexity O(nd log3 ( n )). In Rd , one may hope
that such an algorithm can be converted into an exact
medoid algorithm, but it is not clear how to do this.
Thus, while it may be possible that fast geometric median algorithms can provide inspiration in the development of medoid algorithms, they do not work out
of the box. Moreover, geometric median algorithms
cannot be used for network data as they only work
in vector spaces, thus they are useless for the spatial
network datasets which we consider in §5.
2.2

Medoid Algorithms : TOPRANK and
TOPRANK2

In Eppstein and Wang (2004), the RAND algorithm
for estimating the energy of all elements of a set
S = {x(1), . . . , x(N )} is presented. While RAND is presented in the context of graphs, where the N elements
are nodes of an undirected graph and the metric is
shortest path length, it can equally well be applied to
any set endowed with a distance. The simple idea of
RAND is to estimate the energy of each element from a
sample of anchor nodes I, so that for j ∈ {1, . . . , N },
1 X
Ê(j) =
dist(x(j), x(i)).
|I|
i∈I

James Newling, François Fleuret

An elegant feature of RAND in the context of sparse
graphs is that Dijkstra’s algorithm needs only be run
from anchor nodes i ∈ I, and not from every node.
The key result of Eppstein and Wang (2004) is the
following. Suppose that S has diameter ∆, that is
∆=

max

(i,j)∈{1,...,N }2

dist(x(i), x(j)),

and let  > 0 be some error tolerance. If I is of size
Ω(log(N )/), then P(|E(j) − Ê(j)| > ∆) is O N12
for all j ∈ {1, . . . , N }. Using the union bound, this
means there is a O N1 probability that at least one
energy estimate is off by more than ∆, and so we say
that with high probability (w.h.p.) all errors are less
than ∆.
RAND forms the basis of the TOPRANK algorithm
of Okamoto et al. (2008). Whereas RAND w.h.p. returns an element which has energy within  of the
minimum, TOPRANK is designed to w.h.p. return the
true medoid. In motivating TOPRANK, Okamoto et al.
(2008) observe that the expected difference between
consecutively ranked energies is O(∆/N ), and so if one
wishes to correctly rank all nodes, one needs to distinguish between energies at a scale  = ∆/N , for which
the result of Eppstein and Wang (2004) dictates that
Θ(N log N ) anchor elements are required with RAND,
which is more elements than S contains. However, to
obtain just the highest ranked node should require less
information than obtaining a full ranking of nodes, and
it is to this task that TOPRANK is adapted.
The idea behind TOPRANK is to accurately estimate
only the energies of promising elements. The algorithm proceeds in two passes, where in the first pass
promising elements are earmarked. Specifically, the
first pass runs RAND with N 2/3 log1/3 (N ) anchor elements to obtain Ê(i) for i ∈ {1, . . . , N }, and then
discards elements whose Ê(i) lies below threshold τ
given by,
 31

ˆ 0 log n
τ = arg min Ê(j) + 2∆α
,
(2)
n
j∈{1,...,N }
ˆ is an upper bound on ∆ obtained from the
where ∆
anchor nodes, and α0 is some constant satisfying α0 >
1. The second pass computes the true energy of the
undiscarded elements, returning the one with lowest
true energy. Note that a smaller α0 value results in a
lower (better) threshold, we discuss this point further
in SM-C.
To obtain run time guarantees, TOPRANK requires that
the distribution of node energies is non-decreasing near
to the minimum, denoted by E ∗ . More precisely, letting fE be the probability distribution of energies, the
algorithms require the existence of  > 0 such that,
E ∗ ≤ ẽ < e < E ∗ +  =⇒ fE (ẽ) ≤ fE (e).

(3)

5

If assumption 3 holds, then the run time is Õ(N 3 ). A
second algorithm presented in Okamoto et al. (2008)
is TOPRANK2, where the anchor set I is grown incrementally until some heuristic criterion is met. There
is no runtime guarantee for TOPRANK2, although it has
the potential to run much faster than TOPRANK under
favourable conditions. Pseudocode for RAND, TOPRANK
and TOPRANK2 is presented in SM-C.
2.3

K-medoids algorithm : KMEDS

The Voronoi iteration algorithm, which we refer to as
KMEDS, is similar to lloyd, the main difference being
that cluster medoids are computed instead of cluster
means. It has been desribed in the literature at least
twice, once in Hastie et al. (2001) and then in Park
and Jun (2009), where a novel initialisation scheme is
developed. Pseudocode is presented in SM-B.
All N 2 distances are computed and stored upfront with
KMEDS. Then, at each iteration, KN comparisons are
made during assignment and Ω(N 2 /K) additions are
made during medoid update. The initialisation scheme
of KMEDS requires all N 2 distances. Each iteration of
KMEDS requires retrieving at least max KN, N 2 /K
distinct distances, as can be shown by assuming balanced clusters.
As an alternative to computing all distances upfront,
one could store per-cluster distance matrices which get
updated on-the fly when assignments change. Using
such an approach,the best one could hope for would be
max KN, N 2 /K distance calculations and Θ(N 2 /K)
memory. If one were to completely forego storing distances in memory and calculate distances only when
needed, the number of distance calculations would be
at least r(KN + N 2 /K), where r is the number of
iterations.
The initialisation scheme of Park and Jun (2009) selects K well centered elements as initial medoids. This
goes against the general wisdom for K-means initialisation, where centroids are initialised to be well separated (Arthur and Vassilvitskii, 2007). While the new
scheme of Park and Jun (2009) performs well on a
limited number of small 2-D datasets, we show in § 3
that in general uniform initialisation performs as well
or better.

3

Our new medoid algorithm : trimed

We present our new algorithm, trimed, for determining the medoid of set S = {x(1), . . . , x(N )}. Whereas
the approach with TOPRANK is to empirically estimate
E(i) for i ∈ {1, . . . , N }, the approach with trimed,
presented as Alg. 1, is to bound E(i). When trimed
terminates, an index m∗ ∈ {1, . . . , N } has been de-

A Sub-Quadratic Exact Medoid Algorithm

termined, along with lower bounds l(i) for all i ∈
{1, . . . , N }, such that E(m∗ ) ≤ l(i) ≤ E(i), and thus
x(m∗ ) is the medoid. The bounding approach uses the
triangle inequality, as depicted in Figure 1.
Algorithm 1 The trimed algorithm for computing
the medoid of {x(1), . . . , x(N )}.
1: l ← 0N // lower bounds on energies, maintained
such that l(i) ≤ E(i) and initialised as l(i) = 0.
2: mcl , E cl ← −1, ∞ // index of best medoid candidate found so far, and its energy.
3: for i ∈ shuffle ({1, . . . , N }) do
4:
if l(i) < E cl then
5:
for j ∈ {1, . . . , N } do
6:
d(j) ← dist(x(i), x(j))
7:
end for P
N
8:
l(i) ← N1 j=1 d(j) // set l(i) to be tight,
that is l(i) = E(i).
9:
if l(i) < E cl then
10:
mcl , E cl ← i, l(i)
11:
end if
12:
for j ∈ {1, . . . , N } do
13:
l(j) ← max(l(j), |l(i) − d(j)|) // using
E(i) and dist(x(i), x(j)) to possibly improve bound on E(j).
14:
end for
15:
end if
16: end for
17: m∗ , E ∗ ← mcl , E cl
18: return x(m∗ )
The algorithm trimed iterates through the N elements
of S. Each time a new element with energy lower than
the current lowest energy (E cl ) is found, the index
of the current best medoid (mcl ) is updated (line 10).
Lower bounds on energies are used to quickly eliminate
poor medoid candidates (line 4). Specifically, if lower
bound l(i) on the energy of element i is greater than
or equal to E cl , then i is eliminated. If the bound test
fails to eliminate element i, then it is computed, that is,
all distances to element i are computed (line 6). The
computed distances are used to potentially improve
lower bounds for all elements (line 13). Theorem 3.1
states that trimed finds the medoid. The proof relies
on showing that lower bounds remain consistent when
updated (line 13).
The algorithm is very straightforward to implement,
and requires only two additional floating point values
per datapoint: for sample i, one for l(i) and one for
d(i). Computing either all or no distances from a sample makes particularly good sense for network data,
where computing all distances to a single node is efficiently performed using Dijkstra’s algorithm.
Theorem 3.1. trimed returns the medoid of set S.

x(j)

x(i)

x(j)
x(i)

Figure 1: Using the inequality E(j) ≥ |E(i) −
dist(x(i), x(j)) | to eliminate x(j) as a medoid candidate. Computed element x(i) with energy E(i) ≥ E cl
is used as a pivot to lower bound E(j). The two
cases where the inequality is effective are when (case
1, above) dist(x(i), x(j)) − E(i) ≥ E cl and (case 2,
below) E(i) − dist(x(i), x(j)) ≥ E cl , as both lead to
E(j) ≥ E cl which eliminates x(j) as a medoid candidate.

Proof. We need to prove that l(j) ≤ E(j) for all j ∈
{1, . . . , N } at all iterations of the algorithm. Clearly,
as l(j) = 0 at initialisation, we have l(j) ≤ E(j) at initialisation. E(j) does not change, and the only time
that l(j) may change is on line 13, where we need to
check that |l(i) − d(j)| ≤ E(j). At line 13, l(i) = E(i)
from line 8, and d(j) = dist(x(i), x(j)), so at line 13 we
are effectively checking that |E(i) − dist(x(i), x(j))| ≤
E(j). But this is a simple consequence of the triangle inequality,
PNas we now show. Using the definition,
E(j) = N1 l=1 dist(x(l), x(j)), we have on the one
hand,
E(j) ≥

N
1 X
dist(x(l), x(i)) − dist(x(i), x(j))
N
l=1

≥ E(i) − dist(x(i), x(j)),

(4)

and on the other hand,
E(j) ≥

N
1 X
dist(x(i), x(j)) − dist(x(l), x(i))
N
l=1

≥ dist(x(i), x(j)) − E(i).

(5)

Combining (4) and (5) we obtain the required inequality |E(i) − dist(x(i), x(j))| ≤ E(j).
The bound test (line 4) becomes more effective at later
iterations, for two reasons. Firstly, whenever an element is computed, the lower bounds of other samples
may increase. Secondly, E cl will decrease whenever a
better medoid candidate is found. The main result of
this paper, presented as Theorem 3.2, is that in Rd
1
the expected number of computed elements is O(N 2 )
under some weak assumptions. We show in §5 that

James Newling, François Fleuret
1

The shuffle on line 3 is performed to avoid w.h.p.
pathological orderings, such as when elements are ordered in descending order of energy which would result
in all N elements being computed.
Theorem 3.2. Let S = {x(1), . . . , x(N )} be a set of
N elements in Rd , drawn independently from probability distribution function fX . Let the medoid of S be
x(m∗ ), and let E(m∗ ) = E ∗ . Suppose that there exist strictly positive constants ρ, δ0 and δ1 such that for
any set size N with probability 1 − O(1/N )
x ∈ Bd (x(m∗ ), ρ) =⇒ δ0 ≤ fX (x) ≤ δ1 ,

∗

(7)
∗

2

E(i) − E ≥ αkx(i) − x(m )k .
Then, the expected
number of elements
computed by

  1
4 d
N 2 , where Vd [1] =
trimed is O Vd [1]δ1 + d α
d

π 2 /(Γ( d2 + 1)) is the volume of Bd (0, 1).
3.1

On the assumptions in Theorem 3.2

The assumption of constants ρ, δ0 and δ1 made in Theorem 3.2 is weak, and only pathological distributions
might fail it, as we now discuss. For the assumptions to
fail requires that fX vanishes or diverges at the distribution medoid. Any reasonably behaved distribution
does not have this behaviour, as illustrated in Figure 2.
The constant α is a strong convexity constant. The
existence of α > 0 is guaranteed by the existence of
ρ, δ0 and δ1 , as the mean of a sum of uniformly spaced
cones converges to a quadratic function. This is illustrated in 1-D in Figure 5 in SM-G, but holds true in
any dimension.
Note that the assumptions made are on the distribution fX , and not on the data itself. This must be so
in order to prove complexity results in N .
3.2

fX
δ0

αkx − x(m∗ )k2

(6)

where Bd (x, r) = {x0 ∈ Rd : kx0 − xk ≤ r}. Let
α > 0 be a constant (independent of N ) such that with
probability 1 − O(1/N ) all i ∈ {1, . . . , N } satisfy,
x(i) ∈ Bd (x(m∗ ), ρ) =⇒

δ1

E − E∗

the O(N 2 ) result holds even in settings where the assumptions are not valid or relevent, such as for network
data.

Sketch of proof of Theorem 3.2

We now sketch the proof of Theorem 3.2, showing
how (6) and (7) are used. A full proof is presented
in SM-G. Firstly, let the index of the first element
after the shuffle on line 3 be i0 . Then, no elements
beyond radius 2E(i0 ) of x(i0 ) will subsequently be
computed, due to type 1 eliminations (see Figure 1).
Therefore, all computed elements are contained within
Bd (x(i0 ), 2E(i0 )).

x(m∗ )

x(m∗ ) + ρ

Figure 2: Illustration in 1-D of the constants used in
Theorem 3.2. Above, δ0 and δ1 bound the probability
density function in a region containing the distribution
medoid. Below, the energy of samples grows quadratically around the medoid x(m∗ ). The energy E is a sum
of cones centered on samples, which is approximately
quadratic unless fX vanishes or explodes, guaranteeing the existence of α > 0 required in Theorem 3.2.

Next, notice that once an element x(i) has been
computed in trimed, no elements in the ball
Bd (x(i), E(i) − E cl ) will subsequently be computed,
due to type 2 eliminations (see Figure 1). We refer to
such a ball as an exclusion ball. By upper bounding the
number of exclusion balls contained in Bd (x(i0 ), 2E(i0 ))
using a volumetric argument, we can obtain a bound
on the number of computed elements, but obtaining
such an upper bound requires that the radii of exclusion ball E(i)−E cl be bounded below by a strictly positive value. However, by using a volumetric argument
only beyond a certain positive radius of the medoid
(a radius N −1/2d ), we have α > 0 in (15) which provides a lower bound on exclusion ball radii, assuming
E cl ≈ E ∗ . Using δ0 we can show that E cl approaches
E ∗ sufficiently fastsufficiently fast to validate the approximation E cl ≈ E ∗ .
It then remains to count the number of computed
elements within radius N −1/2d of the medoid. One
cannot find a strict upper bound here, but using the
boundedness of fX provided by δ1 , we have w.h.p. that
the number of elements computed within N −1/2d is
O(δ1 N 1/2 ), as the volume of a sphere scales as the
d’th power of its radius.

A Sub-Quadratic Exact Medoid Algorithm

4

Our accelerated K-medoids
algorithm : trikmeds

We adapt our new medoid algorithm trimed and borrow ideas from Elkan (2003) to show how KMEDS can
be accelerated. We abandon the initial N 2 distance
calculations, and only compute distances when necessary. The accelerated version of lloyd of Elkan
(2003) maintains KN bounds on distances between
points and centroids, allowing a large proportion of
distance calculations to be eliminated. We use this
approach to accelerate assignment in trikmeds, incurring a memory cost O(KN ). By adopting the
algorithm of Newling and Fleuret (2016a) or that
of Hamerly (2010), the memory overhead can be reduced to O(N ). We accelerate the medoid update step
by adapting trimed, reusing lower bounds between iterations, so that trimed is only run from scratch once
at the start. Details and pseudocode are presented in
SM-H.
One can relax the bound test in trimed so that for
 > 0 element i is computed if l(i)(1 + ) < E cl , guaranteeing that an element with energy within a factor
1 +  of E ∗ is found. It is also possible to relax the
bound tests in the assignment step of trikmeds, such
that the distance to an assigned cluster’s medoid is always within a factor 1+ of the distance to the nearest
medoid. We denote by trikmeds- the trikmeds algorithm where the update and assignment steps are
relaxed as just discussed, with trikmeds-0 being exactly trikmeds. The motivation behind such a relaxation is that, at all but the final few iterations, it is
probably a waste of computation obtaining medoids
and assignments at high resolution, as in subsequent
iterations they may change.

5

Results

We first compare the performance of the medoid algorithms TOPRANK, TOPRANK2 and trimed. We then compare the K-medoids algorithms, KMEDS and trikmeds.
5.1

Medoid algorithm results

We compare our new exact medoid algorithm trimed
with state-of-the-art approximate algorithms TOPRANK
and TOPRANK2. Recall, Okamoto et al. (2008) prove
that the approximate algorithms return w.h.p. the true
medoid. We confirm that this is the case in all our experiments, where the approximate algorithms return
the same element as trimed, which we know to be
correct by Theorem 3.1. We now focus on comparing computational costs, which are proportional to the
number of computed points.

Results on artificial datasets are presented in Figure 3,
where our two main observations relate to scaling in
N and dimension d. The artificial data are (left)
uniformly drawn from [0, 1]d and (right) drawn from
Bd (0, 1) with probability of lying within radius 1/21/d
of 1/200, as opposed to 1/2 as would be the case under uniform density. Details about sampling from this
distribution can be found in SM-F. Results on a mix of
publicly available real and artificial datasets are presented in Table 1 and discussed in §5.1.2.
5.1.1

Scaling with N and d on artificial
datasets

In Figure 3 we observe that the number of points computed by trimed is O(N 1/2 ), as predicted by Theorem 3.2. This is illustrated (right) by the close fit of
the number of computed points to exact square root
curves at sufficiently large N for d ∈ {2, 6}.
Recall that TOPRANK consists of two passes, a first
where N 2/3 log1/3 N anchor points are computed, and
a second where all sub-threshold points are computed.
We observe that for small N TOPRANK computes all
N points, which corresponds to all points lying below threshold. At sufficiently large N the threshold
becomes low enough for all points to be eliminated after the first pass. The effect is particularly dramatic
in high dimensions (d = 6 on right), where a phase
transition is observed between all and no points being
computed in the second pass.
Dimension d appears in Theorem 3.2 through a factor
d(4/α)d , where α is the strong convexity of the energy
at the medoid. In Figure 3, we observe that the number of computed points increases with d for fixed N ,
corresponding to a relatively small α. The effect of α
on the number of computed elements is considered in
greater detail in SM-F.
In contrast to the above observation that the number
of computed points increases as dimension increases
for trimed, TOPRANK appears to scale favourably with
dimension. This observation can be explained in terms
of the distribution of energies, with energies close to E ∗
being less common in higher dimensions, as discussed
in SM-J.
5.1.2

Results on publicly available real and
simulated datasets

We present the datasets used here in detail in SM-I.
For all datasets, algorithms TOPRANK, TOPRANK2 and
trimed were run 10 times with a distinct seed, and
the mean number of iterations (n̂) over the 10 runs
was computed. We observe that our algorithm trimed
is the best performing algorithm on all datasets, although in high-dimensions (MNIST-0) and on social

James Newling, François Fleuret

number of computed elements

106

106

TOPRANK
trimed
d=2
d=3
d=4
d=5
d=6

105

104

105

d = 2 TOPRANK
d = 2 trimed
d = 6 TOPRANK
d = 6 trimed

104

103

N 1
18N 2
1
2
N 3 1log 3 N
3N 2

103

102 2
10

103

104
N

105

106

102 3
10

104

105
N

106

107

Figure 3: Comparison of TOPRANK and our algorithm trimed on simulated data. On the left, points are drawn
uniformly from [0, 1]d for d ∈ {2, . . . , 6}, and on the right they are drawn from Bd (0, 1) for d ∈ {2, 6}, with an
increased density near the edge of the ball. Fewer points (elements) are computed by trimed than by TOPRANK
in all scenarios. For small N , TOPRANK computes O(N ) points, before transitioning to Õ(N 2/3 ) computed points
for large N . trimed computes O(N 1/2 ) points. Note that trimed performs better in low-d than in high-d, with
the reverse trend being true for TOPRANK. These observations are discussed in further detail in the text.
network data (Gnutella) no algorithm computes significantly fewer than N elements. The failure in
high-dimensions (MNIST-0) of trimed is in agreement with Theorem 3.2, where dimension appears as
the exponent of a constant term. The small world
network data, Gnutella, can be embedded in a highdimensional Euclidean space, and thus the failure on
this dataset can also be considered as being due to
high-dimensions. For low-dimensional real and spatial
network data, trimed consistently computes O(N 1/2 )
elements.
5.1.3

But who needs the exact medoid
anyway?

A valid criticism that could be raised at this stage
would be that for large datasets, finding the exact
medoid is probably overkill, as any point with energy reasonably close to E ∗ suffices for most applications. But consider, the RAND algorithm requires
computing log N/2 elements to confidently return an
element with energy within E ∗ of E ∗ . For N = 105
and  = 0.05, this is 4600, already more than trimed
requires to obtain the exact medoid on low-d datasets
of comparable size.
5.2

K-medoids algorithm results

With N elements to cluster, KMEDS is Θ(N 2 ) in memory, rendering it unusable on even moderately large
datasets. To compare the initialisation scheme proposed in Park and Jun (2009) to random initialisation,
we have performed experiments on 14 small datasets,
with K ∈ {10, dN 1/2 e, dN/10e}. For each of these 42

experimental set-ups, we run the deterministic KMEDS
initialisation once, and then uniform random initialisation, 10 times. Comparing the mean final energy
of the two initialisation schemes, in only 9 of 42 cases
does KMEDS initialisation result in a lower mean final
energy. A Table containing all results from these experiments in presented in SM-E.
Having demonstrated that random uniform initialisation performs at least as well as the initialisation
scheme of KMEDS, and noting that trikmeds-0 returns
exactly the same clustering as would KMEDS with uniform random initialisation, we turn our attention to
the computational performance of trikmeds. Table 2
presents results on 4 datasets, each described in SMI. The first numerical column is the relative number
of distance calculations using trikmeds-0 and KMEDS,
where large savings in distance calculations, especially
in low-dimensions, are observed. Columns φc and φE
are the number of distance calculations and energies
respectively, using  ∈ {0.01, 0.1}, relative to  = 0.
We observe large reductions in the number of distance
computations with only minor increases in energy.

6

Conclusion and future work

We have presented our new trimed algorithm for computing the medoid of a set, and provided strong theoretical guarantees about its performance in Rd . In
low-dimensions, it outperforms the state-of-the-art approximate algorithm on a large selection of datasets.
The algorithm is very simple to implement, and can
easily be extended to the general ranking problem.
In the future, we propose to explore the idea of us-

A Sub-Quadratic Exact Medoid Algorithm

dataset
Birch 1
Birch 2
Europe
U-Sensor Net
D-Sensor Net
Pennsylvania road
Europe rail
Gnutella
MNIST

type
2-d
2-d
2-d
u-graph
d-graph
u-graph
u-graph
d-graph
784-d

N
1.0 × 105
1.0 × 105
1.6 × 105
3.6 × 105
3.6 × 105
1.1 × 106
4.6 × 104
6.3 × 103
6.7 × 103

TOPRANK TOPRANK2 trimed
n̂
n̂
n̂
57944
100180
2180
66062
100180
2208
176095
169535
2862
113838
327216
1593
99896
176967
1372
216390
time-out
2633
35913
47041
518
7043
6407
6328
7472
6799
6514

Table 1: Comparison of TOPRANK, TOPRANK2 and our algorithm trimed on publicly available real and simulated
datasets. Column 2 provides the type of the dataset, where ‘x-d’ denotes x-dimensional vector data, while
‘d-graph’ and ‘u-graph’ denote directed and undirected graphs respectively. Column n̂ gives the mean number
of elements computed over 10 runs. Our proposed trimed algorithm obtains the true medoid with far fewer
computed points in low dimensions and on spatial network data. On the social network dataset (Gnutella) and
the very high-d dataset (MNIST), all algorithms fail to provide speed-up, computing approximately N elements.

Dataset
Europe
Conflong
Colormo
MNIST50

N
1.6 × 105
1.6 × 105
6.8 × 104
6.0 × 104

=0
d Nc /N 2
2 0.067
3 0.042
9 0.163
50 0.280

K = 10
 = 0.01
φc
φE
0.33 1.004
0.67 1.001
0.92 1.000
0.99 1.000

 = 0.1
=0
φc
φE Nc /N 2
0.01 1.054 0.008
0.08 1.014 0.006
0.35 1.015 0.011
0.95 1.001 0.019

√
K = d Ne
 = 0.01
 = 0.1
φc
φE
φc
φE
0.68 1.031 0.39 1.090
0.92 1.003 0.61 1.026
0.98 1.000 0.82 1.005
0.99 1.001 0.97 1.001

Table 2: Relative numbers of distance calculations and final energies using trikmeds- for  ∈ {0, 0.01, 0.1}. The
number of distance calculations with trikmeds-0 is Nc , presented here relative to the number computed using
KMEDS (N 2 ) in column Nc /N 2 . The number of distance calculations with  ∈ {0.01, 0.1} relative to trikmeds-0
are given in columns φc , so φc = 0.33 means 3× fewer calculations than with  = 0. The final energies with
 ∈ {0.01, 0.1} relative to trikmeds-0 are given in columns φE . We see that trikmeds-0 uses significantly fewer
distance calculations than would KMEDS, especially in low-dimensions where a greater than K× reduction is
observed (NC /N 2 < 1/K). For low-d, additional relaxation further increases the saving in distance calculations
with little cost to final energy.

ing more complex triangle inequality bounds involving
several points, with as goal to improve on the O(N 1/2 )
number of computed points.
We have demonstrated how trimed, when combined
with the approach of Elkan (2003), can greatly reduce
the number of distance calculations required by the
Voronoi iteration K-medoids algorithm of Park and
Jun (2009). In the future we would like to replace the
strategy of Elkan (2003) with that of Hamerly (2010),
which will be better adapted to graph clustering as
either all or no distances are computed with it, making
it more amenable to Dijkstra’s algorithm.

Acknowledgements
The authors are grateful to Wei Chen for helpful discussions of the TOPRANK algorithm. James Newling
was funded by the Hasler Foundation under the grant
13018 MASH2.

References
Arthur, D. and Vassilvitskii, S. (2007). K-means++:
The advantages of careful seeding. In Proceedings of
the Eighteenth Annual ACM-SIAM Symposium on
Discrete Algorithms, SODA ’07, pages 1027–1035,
Philadelphia, PA, USA. Society for Industrial and
Applied Mathematics.
Chipman, H., Hastie, T., and Tibshirani, R. (2003).

James Newling, François Fleuret

Statistical Analysis of Gene Expression Microarray
Data. Chapman & Hall. Chapter 4.
Cohen, E., Delling, D., Pajor, T., and Werneck, R. F.
(2014). Computing classic closeness centrality, at
scale. In Proceedings of the Second ACM Conference
on Online Social Networks, COSN ’14, pages 37–50,
New York, NY, USA. ACM.
Cohen, M. B., Lee, Y. T., Miller, G. L., Pachocki,
J. W., and Sidford, A. (2016). Geometric median in
nearly linear time. In STOC16. submitted.
Elkan, C. (2003). Using the triangle inequality to accelerate k-means. In Machine Learning, Proceedings
of the Twentieth International Conference (ICML
2003), August 21-24, 2003, Washington, DC, USA,
pages 147–153.
Eppstein, D. and Wang, J. (2004). Fast approximation
of centrality. J. Graph Algorithms Appl., 8(1):39–45.
Frahm, J.-M., Fite-Georgel, P., Gallup, D., Johnson,
T., Raguram, R., Wu, C., Jen, Y.-H., Dunn, E.,
Clipp, B., Lazebnik, S., and Pollefeys, M. (2010).
Building rome on a cloudless day. In Proceedings of
the 11th European Conference on Computer Vision:
Part IV, ECCV’10, pages 368–381, Berlin, Heidelberg. Springer-Verlag.
Hamerly, G. (2010). Making k-means even faster. In
SDM, pages 130–140.
Hastie, T. J., Tibshirani, R. J., and Friedman, J. H.
(2001). The elements of statistical learning : data
mining, inference, and prediction. Springer series in
statistics. Springer, New York.
Hoare, C. A. R. (1961). Algorithm 65: Find. Commun.
ACM, 4(7):321–322.
Kaufman, L. and Rousseeuw, P. J. (1990). Finding
groups in data : an introduction to cluster analysis. Wiley series in probability and mathematical
statistics. Wiley, New York. A Wiley-Interscience
publication.
Newling, J. and Fleuret, F. (2016a). Fast k-means
with accurate bounds. In Proceedings of the International Conference on Machine Learning (ICML),
pages 936–944.
Newling, J. and Fleuret, F. (2016b). K-medoids for
k-means seeding. arXiv:1609.04723. Under review.
Ng, R. T., Han, J., and Society, I. C. (2005). Clarans:
A method for clustering objects for spatial data mining. IEEE Transactions on Knowledge and Data
Engineering, pages 1003–1017.

Okamoto, K., Chen, W., and Li, X.-Y. (2008). Ranking of closeness centrality for large-scale social networks. In Proceedings of the 2Nd Annual International Workshop on Frontiers in Algorithmics, FAW
’08, pages 186–195, Berlin, Heidelberg. SpringerVerlag.
Park, H.-S. and Jun, C.-H. (2009). A simple and
fast algorithm for k-medoids clustering. Expert Syst.
Appl., 36(2):3336–3341.
Rattigan, M. J., Maier, M., and Jensen, D. (2007).
Graph clustering with network structure indices. In
Proceedings of the 24th International Conference on
Machine Learning, ICML ’07, pages 783–790, New
York, NY, USA. ACM.

A Sub-Quadratic Exact Medoid Algorithm

SM-A

On the difficulty of the medoid problem

We construct an example showing that no general purpose algorithm exists to solve the medoid problem in
o(N 2 ). Consider an almost fully connected graph containing N = 2m + 1 nodes, where the graph is exactly m
edges short of being fully connected: one node has 2m edges and the others have 2m − 1 edges. The graph has
2m2 edges. With the shortest path metric, it is easy to see that the node with 2m edges is the medoid, hence the
medoid problem is as difficult as finding the node with 2m edges. But, supposing that the edges are provided as
an unsorted adjacency list, it is clearly an O(m2 ) task to determine which node has 2m edges as one must look
at all edges until a node with 2m edges is found. Thus determining the medoid is O(m2 ) which is O(N 2 ).

SM-B

KMEDS pseudocode

Alg. 2 presents the KMEDS algorithm of Park and Jun (2009), with the novel initialisation of KMEDS on line 1.
KMEDS is essentially lloyd, with medoids instead of means.
Algorithm 2 KMEDS for clustering data {x(1), . . . , x(N )} around K medoids
P
1: Set all distances D(i, j) ← kx(i) − x(j)k and sums S(i) ← j∈{1,...,N } D(i, j)
P
2: Initialise medoid indices as K indices minimising f (i) = j∈{1,...,N } D(i, j)/S(j)
3: while Some convergence criterion has not been met do
4:
Assign each element to the cluster whose medoid is nearest to the element
5:
Update cluster medoids according to assignments made above
6: end while

SM-C

RAND, TOPRANK and TOPRANK2 pseudocode

We present pseudocode for the RAND, TOPRANK and TOPRANK2 algorithms of Okamoto et al. (2008), and discuss
the explicit and implicit constants.
SM-C.1

2

1

On the number of anchor elements in TOPRANK : the constant in Θ(N 3 (log N ) 3 )

Note that the number of anchor points used in TOPRANK does not affect the result that the medoid is w.h.p.
1
returned. However, Okamoto et al. (2008) show that by choosing the size of the anchor set to be q (log N ) 3 for
any q, the run time is guaranteed to be Õ(N 5/3 ). They do not suggest a specific q, the optimal q being dataset
dependant. We choose q = 1.
Consider Figure 3 in Section 5.1 for example, where q = 1. Had q be chosen to be less than 1, the line
ncomputed = N 2/3 log1/3 N to which TOPRANK runs parallel for large N would be shifted up or down by log q,
however the N at which the transition from ncomputed = N 2/3 log1/3 N to ncomputed = N 2/3 log1/3 N takes
place would also change.
SM-C.2

On the parameter α0 in TOPRANK and TOPRANK2

The threshold τ in (2) is proportional to the parameter α0 . In Okamoto et al. (2008), it is stated that α0 should
be some value greater than 1. Note that the smaller α0 is, the lower the threshold is, and hence fewer the number
of computed points is, thus α0 = 1.00001 would be a fair choice. We use α0 = 1 in our experiments, and observe
that the correct medoid is returned in all experiments.
Personal correspondence with the authors of Okamoto et al. (2008) has brought into doubt the proof of the
result that the medoid is w.h.p. returned for any α0 where α0 > 1. In our most recent correspondence,
the
√
authors suggest that the w.h.p. result can be proven with the more conservative bound of α0 > 1.5. Moreover,
0
we show in SM-D that α0 > 1 is good enough to return the medoid with probability N −(α −1) , a probability
which still tends to 0 as N grows large, but not a w.h.p. result. Please refer to SM-D for further details on our
correspondence with the authors.

James Newling, François Fleuret

SM-C.3

On the parameters specific to TOPRANK2

In addition to α0 , TOPRANK2 requires two parameters to be set. The first is l0 , the starting anchor set size, and
the second is q, the amount by which l should be incremented at each iteration. Okamoto et al. (2008) suggest
taking l0 to be the number of top ranked nodes required, which in our case would be l0 = k = 1. However, in
our experience this is too small as all nodes lie well within the threshold and thus when l increases there is no
change to number below threshold, which makes the algorithm break out of the search for the optimal l too early.
Indeed, l0 needs to be chosen so that at least some points
√ have energies greater than the threshold, which in our
experiments is already quite large. We choose l0 = N , as any value larger than N 2/3 would make TOPRANK2
redundant to TOPRANK. The parameter q we take to be log N as suggested by Okamoto et al. (2008).
Algorithm 3 RAND for estimating energies of elements of set S (Eppstein and Wang, 2004).
I ← random uniform sample from {1, . . . , N }
// Compute all distances from anchor elements (I), using Dijkstra’s algorithm on graphs
for i ∈ I do
for j ∈ {1, . . . , N } do
d(i, j) ← kx(i) − x(j)k,
end for
end for
// Estimate energies as mean distances to anchor elements
for j ∈ {1, . . .P
, N } do
1
Ê(j) ← |I|
i∈I d(i, j)
end for
return Ê
Algorithm 4 TOPRANK for obtaining top k ranked elements of S (Okamoto et al., 2008).
2

1

1

l ← N 3 (log N ) 3 // Okamoto et al. (2008) state that l should be Θ((log N ) 3 ), the choice of 1 as the constant
is arbitrary (see comments in the text of Section SM-C.1).
Run RAND with uniform random I of size l to get Ê(i) for i ∈ {1, . . . , N }.
Sort Ê so that Ê[1] ≤ Ê[2] ≤ . . . ≤ Ê[N ]
ˆ ← 2 mini∈I maxj∈{1,...,N } kx(i) − x(j)k // where kx(i) − x(j)k computed in RAND
∆


q
Q ← i ∈ {1, . . . , N } | Ê(i) ≤ Ê[k] + 2α0 ∆ log(n)
.
l
Compute exact energies of all elements in Q and return the element with the lowest energy.

SM-D

On the proof that TOPRANK returns the medoid with high probability

Through correspondence with the authors of Okamoto et al. (2008), we have located a small problem in the
proof that the medoid is returned w.h.p. for α0 > 1, the problem lying in the second inequality of Lemma 1. To
arrive at this inequality, the authors have used the fact that for all i,
P(E(i) ≥ Ê(i) + f (l) · ∆) ≥ 1 −

1
,
2N 2

(8)

which is a simple consequence of the Hoeffding inequality as shown in Eppstein and Wang (2004). Essentially (8)
says that, for a fixed node i, from which the mean distance to other nodes is E(i), if one uniformly samples l
distances to i and computes the mean Ê(i), the probability that Ê(i) is less than E(i) + f (l) is greater than
1 − 2N1 2 .
The inequality (8) is true for a fixed node i. However, it no longer holds if i is selected to be the node with the
lowest Ê(i). To illustrate this, suppose that E(i) = 1 for all i, and compute Ê(i) for all i. Let Ê ∗ = arg mini Ê(i).
Now, we have a strong prior on Ê ∗ being significantly less than 1, and (8) no longer holds as a statement made
about Ê ∗ .
In personal correspondence, the authors show that the problem can be fixed by the use of an additional layer
of union bounding, with a correction to be published (if not already done so at time of writing). However, the

A Sub-Quadratic Exact Medoid Algorithm

Algorithm 5 TOPRANK2 for obtaining top k ranked elements of S (Okamoto et al., 2008).
// In Okamoto et al. (2008), it is suggested that l0 be taken as k, which in the case of the medoid problem is
1. We have experimented with several choices for l0 , as discussed in the text.
l ← l0
Run RAND with uniform random I of size l to get Ê(i) for i ∈ {1, . . . , N }.
ˆ ← 2 mini∈I maxj∈{1,...,N } kx(i) − x(j)k // where kx(i) − x(j)k computed in RAND
∆
Sort Ê

 so that Ê[1] ≤ Ê[2] ≤ . . . ≤ Ê[N ] q
log(n)
0
.
Q ← i ∈ {1, . . . , N } | Ê(i) ≤ Ê[k] + 2α ∆
l
g←1
while g is 1 do
p ← |Q|
// The recommendation for q in Okamoto et al. (2008) is log(n), we follow the suggestion
Increment I with q new anchor points
Update Ê for all data according to new anchor points
l ← |I|
ˆ ← 2 mini∈I maxj∈{1,...,N } kx(i) − x(j)k
∆
Sort Ê

 so that Ê[1] ≤ Ê[2] ≤ . . . ≤ Ê[N ] q
log(n)
0
Q ← i ∈ {1, . . . , N } | Ê(i) ≤ Ê[k] + 2α ∆
l
p0 ← |Q|
if p − p0 < log (n) then
g←0
end if
end while
Compute exact energies of all elements in Q and return the element with the lowest energy
additional layer of union bound requires a more conservative constraint
on α0 , which is α0 > 2, although the
√
0
authors propose that the w.h.p. result can be proven
√ with α > 1.5 for N sufficiently large. We now present
a small proof proving the w.h.p. result for α0 > 2 for N sufficiently large, with at the same time α0 > 1
0
guaranteeing that the medoid is returned with probability O(N α −1 ).
SM-D.1

That the medoid is returned with high probability holds for α0 >
vanishing probability it is returned for α0 > 1

√

2 and that with

Recall that we have N nodes with energies E(1), . . . , E(n). We wish to find the k lowest energy nodes (the
original setting of Okamoto et al. (2008)). From Hoeffding’s inequality we have,

P(|E(i) − Ê(i)| ≥ ∆) ≤ 2 exp −l2 .
(9)
Set the probability on the right hand side of 9 to be 2/N 1+β , that is,

2 exp −l2 = 2/N 1+β ,
which corresponds to
s
=
Clearly

√

1+β
l



log (N ) := f˜(l).

1 + β corresponds to α0 . With this notation we have,
P(|E(i) − Ê(i)| ≥ f˜(l)∆) ≤

2
N 1+β

.

Applying the union bound to (10) we have,
 

2
P ¬ ∧i∈{1,...,N } |E(i) − Ê(i)| ≤ f˜(l)∆ ≤ β .
N

(10)

(11)

James Newling, François Fleuret

K = 10
Dataset
gassensor
house16H
S1
S2
S3
S4
A1
A2
A3
thyroid
yeast
wine
breast
spiral

K=

l√ m
N

K=

N 
10

N
d µu /µpark σu /µpark µu /µpark σu /µpark µu /µpark σu /µpark
256 128
1.09
0.08
0.90
0.03
0.83
0.01
1927 17
1.01
0.02
0.97
0.01
0.93
0.01
5000 2
1.05
0.05
0.75
0.01
0.32
0.01
5000 2
1.04
0.07
0.68
0.01
0.34
0.00
5000 2
1.03
0.05
0.76
0.01
0.35
0.00
5000 2
1.02
0.03
0.75
0.01
0.41
0.01
3000 2
0.82
0.03
0.43
0.01
0.19
0.00
5250 2
0.98
0.03
0.47
0.01
0.25
0.00
7500 2
0.96
0.02
0.42
0.02
0.22
0.00
215 5
0.95
0.08
0.97
0.04
0.93
0.04
1484 8
1.00
0.02
0.96
0.02
0.91
0.02
178 14
1.01
0.02
1.02
0.01
0.98
0.02
699 9
0.79
0.03
0.77
0.02
0.68
0.02
312 3
1.03
0.03
0.99
0.02
0.82
0.03

Table 3: Comparing the initialisation scheme proposed in Park and Jun (2009) with random uniform initialisation
for the KMEDS algorithm. The final energy using the deterministic scheme proposed in Park and Jun (2009) is
µpark . The mean over 10 random uniform initialisations is µu , and the corresponding standard deviation is σu .
For small K (K = 10), the performances using the two schemes are comparable, while for larger K, it is clear
that uniform initialisation performs much better on the majority of datasets.
Recall that we wish to obtain the k nodes with lowest energy. Denote by r(j) the index of the node with the
j’th lowest energy, so that
E(r(1)) ≤ . . . ≤ E(r(j)) ≤ . . . ≤ E(r(N )).
Denote by r̂(j) the index of the node with the j’th lowest estimated energy, so that
Ê(r̂(1)) ≤ . . . ≤ Ê(r̂(j)) ≤ . . . ≤ Ê(r̂(N )).
Now assume that for all i, it is true that |E(i) − Ê(i)| ≤ f˜(l). Then consider, for j ≤ k,
 


 
Ê(r̂(k)) − Ê(r(j)) = Ê(r̂(k)) − E(r(k)) + E(r(k)) − E(r(j)) + E(r(j)) − Ê(r(j)) ,
|
{z
} |
{z
} |
{z
}
≥−f˜(l)∆

≥0

(12)

≥−f˜(l)∆

≥ −2f˜(l)∆.
The first bound in (12) is obtained by considering the most extreme case possible under the assumption, which
is Ê(i) = a(E) − f˜(l) for all i. The second bound follows from j ≤ k, and the third bound follows directly from
the assumption. We thus have that, under the assumption,
Ê(r(j)) ≤ Ê(r̂(k)) + 2f˜(l)∆,
which says that all nodes of rank less than or equal to k have approximate energy less than Ê(r̂(k)) + 2f˜(l)∆.
As the assumption holds with probability greater than 1 − 2/N β by (11), we are done. Take β = 1 if you want
the statement with high probability, that is
r
2 log(n)
,
=
l
but for any β > 0, which corresponds to α0 > 1, the probability of failing to return the k lowest energy nodes
tends to 0 as N grows.

SM-E

On the initialisation of Park and Jun (2009)

In Table 3 we present the full results of the 48 experiments comparing the initialisation proposed in Park and
Jun (2009) with simple uniform initialisation. The 14 datasets are all available from https://cs.joensuu.fi/
sipu/datasets/.

number of points computed

A Sub-Quadratic Exact Medoid Algorithm

2.5 1e4
2.0
1.5

d =2
d =3
d =4
d =5

p

22.3 N
p

15.5 N

1.0

p

10.6p N
7.8 N

0.5
0.0
0.0

0.2

0.4

0.6
N

0.8

1.0

1.2 1e4
1.0
0.8
0.6
0.4
0.2
0.0
1.2 0.0
0.2
1e6

p

11.2 N
p

7.1 N
p

5.0p N
3.7 N

0.4

0.6
N

0.8

1.0

1.2
1e6

Figure 4: Number of points computed on simulated data. Points are drawn from Bd (0, 1), for d ∈ {2, 3, 4, 5}.
On the left, points are drawn uniformly, while on the right, the density in Bd (0, (1/2)1/d ) is 19× lower that
in Ad (0, (1/2)1/2 , 1), where recall that Ad (x, r1 , r2 ) denotes an annulus centred at x√of inner radius r1 and
outer radius r2 . We observe a near perfect fit of the number of computed points to ξ N where the constant
ξ depends on the dimension and the distribution (left and right). The number of computed points increases
with dimension. The strong convexity constant of the distribution on the right is larger, corresponding to fewer
distance calculations as predicted by Theorem 3.2.

SM-F

Scaling with α, N , and dimension d

We perform more experiments to provide further validation of Theorem 3.2. In particular, we check how the
number of computed elements scales with N , d and α. We generate data from a unit ball in various dimensions,
according to two density functions with different strong convexity constants α. The first density function is
uniform, so that the density everywhere in the ball is uniform. To sample from this distribution, we generate
two random variables, X1 ∼ Nd (0, 1) and X2 ∼ U (0, 1) and use
1

X3 = X1 /kX1 k · X2d ,

(13)

as a sample from the unit ball Bd (0, 1) with uniform distribution. The second distribution we consider has a
higher density beyond radius (1/2)1/d . Specifically, within this radius the density is 19× lower than beyond this
radius. To sample from this distribution, we sample X3 according to (13), and then points lying within radius
(1/2)1/d are with probability 1/10 re-sampled uniformly beyond this radius.
The second distribution has a larger strong convexity constant α. To see this, note that the strong convexity
constant at the center of the ball depends only on the density of the ball on its surface, that is at radius 1, as can
be shown using an argument based on cancelling energies of internal points. As the density at the surface under
distribution 2 is approximately twice that of under distribution 1, the change in energy caused by a small shift
in the medoid is twice as large under distribution 2. Thus, according to Theorem 3.2, we expect the number of
computed points to be larger under distribution 1 than under distribution 2. This is what we observe, as shown
in Figure 4, where distribution 1 is on the left and distribution 2 is on the right.
In Figure 4 we observe a near perfect N 1/2 scaling of number of computed points. Dashed curves are exact N 1/2
relationships, while the coloured points are the observed number of computed points.

SM-G

Proof of Theorem 3.2 (See page 5)

Theorem 3.2. Let S = {x(1), . . . , x(N )} be a set of N elements in Rd , drawn independently from probability
distribution function fX . Let the medoid of S be x(m∗ ), and let E(m∗ ) = E ∗ . Suppose that there exist strictly
positive constants ρ, δ0 and δ1 such that for any set size N with probability 1 − O(1/N )
x ∈ Bd (x(m∗ ), ρ) =⇒ δ0 ≤ fX (x) ≤ δ1 ,

(6)

James Newling, François Fleuret

where Bd (x, r) = {x0 ∈ Rd : kx0 −xk ≤ r}. Let α > 0 be a constant (independent of N ) such that with probability
1 − O(1/N ) all i ∈ {1, . . . , N } satisfy,
x(i) ∈ Bd (x(m∗ ), ρ) =⇒

(7)

∗

∗

2

E(i) − E ≥ αkx(i) − x(m )k .

Then, the expected number of elements computed by trimed is O Vd [1]δ1 + d


4 d
α




1
N 2 , where Vd [1] =

d

π 2 /(Γ( d2 + 1)) is the volume of Bd (0, 1).
Proof. We show that the assumptions made in Th. 3.2 validate the assumptions required in Thm SM-G.1.
Firstly, if e(i) > ρ then e(i) ≥ αρ2 e(i) > ρ, which follows from the convexity of the loss function and. Secondly,
the existance of β follows from continuity of the gradient of the distance, combined with the existence of δ1
(non-exploding).
Theorem SM-G.1 (Main Theorem Expanded). Let S = {x(1), . . . , x(N )} ⊂ Rd have medoid x(m∗ ) with minimum energy E(m∗ ) = E ∗ , where elements in S are drawn independently from probability distribution function
fX . Let e(i) = kx(i) − x(m∗ )k. Suppose that for fX there exist strictly positive constants α, β, ρ, δ0 and δ1
satisfying,
x ∈ Bd (x(m∗ ), ρ) =⇒ δ0 ≤ fX (x) ≤ δ1 ,
(14)
where Bd (x, r) = {x0 ∈ Rd : kx0 − xk ≤ r}, and that for any set size N , w.h.p. all i ∈ {1, . . . , N } satisfy,
(
αe(i)2 if e(i) ≤ ρ,
E(i) − E ∗ ≥
αρ2
if e(i) > ρ,

(15)

and,
E(i) − E ∗ ≤ βe(i)2

if

e(i) ≤ ρ.

(16)

Then
which is to say not eliminated on line 4 of trimed, is
 the expected number
 1  of elements computed,
d
d
O Vd [1]δ1 + d α4
N 2 , where Vd [1] = π 2 /(Γ( d2 + 1)) is the volume of Bd (0, 1).
E − E∗

x(m∗ )

Figure 5: A sum of uniformly distributed cones is approximately quadratic.
1

1

Proof. We first show that the expected number of computed elements in Bd (x(m∗ ), N − 2d ) is O(Vd [1]δ1 N 2 ).
1
When N is sufficiently large, fX (x) ≤ δ1 within Bd (x(m∗ ), N − 2d ). The expected number of samples in
1
Bd (x(m∗ ), N − 2d ) is thus upper bounded by δ1 multiplied by the volume of the ball. But the volume of a
1
1
ball of radius N − 2d in Rd is Vd [1]N − 2 .
In Lemma SM-G.2 we usea packing argument
to show that the number of computed elements in the annulus


1
1
4 d
∗
− 2d
Ad (x(m ), N
, ∞) is O d α N 2 , but we there assume that the medoid index m∗ is the first element in
shuffle({1, . . . , N }) on line 3 of trimed and thus that the medoid energy is known from the first iteration
(E cl = E ∗ ). We now extend Lemma SM-G.2 to the case where the medoid is not the first element processed. We
1
do this by showing that w.h.p. an element with energy very close to E ∗ has been computed after N − 2 iterations
of trimed, and thus that the bounds on numbers of computed elements obtained using the packing arguments
1
underlying Lemma SM-G.2 are all correct to within some small factor after N − 2 iterations.

A Sub-Quadratic Exact Medoid Algorithm

E(i) − E ∗

0.12

0.06

αρ2
0.4

−ρ

ρ

0.0
kx(i) − x(m∗ )k

0.4

Figure 6: Illustrating the parameters α, β and ρ of Theorem 3.2. Here we draw N = 101 samples uniformly from
[−1, 1] and compute their energies, plotted here as the series of points. Theorem 3.2 states that their exists α,
β and ρ such that irrespective of N , all energies (points) will lie in the envelope (non-hatched region).

2

2

The probability of a sample lying within radius N − 3d of x(m∗ ) is Ω(δ0 N − 3 ), and so the probability that none
1

2

2

1
2

1

of the first N 2 samples lies within radius N − 3d is O((1 − δ0 N − 3d )N ) which is O( N1 ). Thus w.h.p. after N 2
4
iterations of trimed, E cl is within βN − 3d of E ∗ , which means that the radii of the balls used in the packing
1
argument are overestimated by at most a factor N − 3d . Thus w.h.p. the upper bounds obtained with the packing
1
argument are correct to within a factor 1 + N − 3 . The remaining O( N1 ) cases do not affect the expectation, as
we know that no more than N elements can be computed.

Lemma SM-G.2 (Packing beyond the vanishing radius). If we assume (15) from Theorem 3.2 and that
the medoid index m∗ is the
element
processed by trimed, then the number of elements computed in
 first


1
1
4 d
∗
− 2d
2
Ad (x(m ), N
, ∞) is O d α N .

Proof. Follows from Lemmas SM-G.3 and SM-G.4.

1

Lemma SM-G.3 (Packing from the vanishing radius N − d to ρ). If we assume (15) from Theorem 3.2 and
that the medoid index m∗ is the first element processed in trimed, then the number of computed elements in
d 1
1
A(x(m∗ ), N − 2d , ρ) is O(d α4 N 2 ).

Proof. According to Assumption 15, an element at radius r < ρ has surplus energy at least αr2 . This means
that, assuming that the medoid has already been computed, an element computed at radius r will be surrounded
by an exclusion zone of radius αr2 in which no element will subsequently be computed. We will use this fact
1
to upper bound the number of computed elements in A(x(m∗ ), N − 2d , ρ), firstly by bounding the number in an
annulus of inner radius r and width αr2 , that is the annulus Ad (x(m∗ ), r, r + αr2 ), and then summing over
1
concentric rings of this form which cover A(x(m∗ ), N − 2d , ρ). Recall that the number of computed elements in
Ad (x(m∗ ), r, r + αr2 ) is denoted by Nc (x(m∗ ), r, r + αr2 ).

James Newling, François Fleuret

We use Lemma SM-G.5 to bound Nc (x(m∗ ), r, r + αr2 ),

d 2
4
αr (r + αr2 )d−1
∗
2
2
Nc (x(m ), r, r + αr ) ≤ (d + 1) √
d
3
(αr2 )

d 
d−1
4
1
≤ (d + 1)2 √
1+
αr
3

d 

d−1
4
2
2
√
≤ (d + 1)
max 2,
αr
3
 d−1 !!
d

2
4
d−1
2
max 2 ,
≤ (d + 1) √
αr
3

d
 d−1 !
4
2
≤ (d + 1)2 √
2d−1 +
αr
3
d

d  d−1

8
1
8
+ (d + 1)2 √
≤ (d + 1)2 √
αr
3
3
1

Let r0 = N − 2d and ri+1 = ri + αri2 , and let T be the smallest index i such that ri ≤ ρ. With this notation in
hand, we have
T
X
1
Nc (x(m∗ ), N − 2d , ρ) ≤
Nc (x(m∗ ), ri , αri + ri2 ).
i=0

The summation on the right-hand side can be upper-bounded by an integral. Using that the difference between
ri and ri+1 is αri2 , we need to divide terms in the sum by αri2 when converting to an integral. Doing this, we
obtain,
Z ρ+αρ2
1
Nc (x(m∗ ), N − 2d , ρ) ≤
Nc (x(m∗ ), r, αr2 )dr
1
N − 2d

d  d Z ∞
8
1
≤ const + (d + 1) √
r−(1+d) dr
1
α
3
N − 2d
 d
1
4
≤ const + (d + 1)
N 2.
α
d
This completes the proof, and provides the hidden constant of complexity as (d + 1) α4 . Thus larger values
for α should result in fewer computed elements in the annulus Ad (x(m∗ ), r, r + αr2 ), which makes sense given
that large values of α imply larger surplus energies and thus larger elimination zones.
2



Lemma SM-G.4 (Packing beyond ρ). If we assume (15) from Theorem 3.2 and that the medoid index m∗ is
the first element processed by trimed, then the number of computed elements in Ad (x(m∗ ), ρ, ∞) is less than
(1 + 4E ∗ /(αρ2 ))d .
Proof. Recall that we at assuming m∗ = 1, that is that the medoid is the first element processed in trimed. All
elements beyond radius 2E ∗ are eliminated by type 1 eliminations (Figure 1), which provides the first inequality
below. Then, as the excess energy is at least  = αρ2 for all elements beyond radius ρ of x(m∗ ), we apply
Lemma SM-G.8 with  = αρ2 /2 to obtain the second inequality below,
Nc (m(x), ρ, ∞) ≤ Nc (m(x), ρ, 2E ∗ )
(2E ∗ + 21 αρ2 )d
( 12 αρ2 )d

d
4E ∗
≤ 1+
.
αρ2
≤

A Sub-Quadratic Exact Medoid Algorithm

Lemma SM-G.5 (Annulus packing). For 0 ≤ r and 0 <  ≤ w. If
X ⊂ Ad (0, r, r + w),
where
∀x ∈ X , Bd (x, ) ∪ X = {x},

(17)

then,
|X | ≤ (d + 1)

2



4
√
3

d

d−1

w (r + w)
d

.

Proof. The condition (17) implies,
 
 
∀x, x0 ∈ X × X , B x,
∪ B x0 ,
= ∅.
2
2

(18)

Using that  ∈ (0, w] and Lemma SM-G.6, one can show that for all x ∈ A(0, r, r + w),
  d2
  

hi
3
1
volume B x,
∩ A(0, r, r + w) >
Vd
2
d+1 4
2

(19)

Combining (18) with (19) we have,
[

volume

x∈X

√ !d
3
|X |d .
4

!

Vd [1]
B x,
∩ A(0, r, r + w) >
2
d+1


(20)

Letting Sd [] denote the surface area of a B(0, ), it is easy to see that
volume (A(0, r, r + w)) < Sd [1] w (r + w)

d−1

.

Combining (20) with (21) we get,
Vd [1]
d+1

√ !d
3
d−1
|X |d < Sd [1] w (r + w)
.
4

which combined with the fact that
Sd [1]
=
Vd [1]

dVd
dr

Vd

!
r=1

= d,
provides us with,
|X | ≤ (d + 1)

2



4
√
3

d

d−1

w (r + w)
d

.

Lemma SM-G.6 (Volume of ball intersection). For x0 , x1 ∈ Rd with kx0 − x1 k = 1,
volume (Bd (x0 , 1) ∩ Bd (x1 , 1))
1
≥
volume (Bd (x0 , 1))
d+1

  d2
3
.
4

(21)

James Newling, François Fleuret

Proof. Let Vd [r] denote the volume of Bd (0, r). It is easy to see that,
Z

1
2

i
hp
x(2 − x) dx
0
"r #
Z 21
3
≥2
Vd−1
x dx
2
0
Z 21   d−1
2
3
dx
≥ 2Vd−1 [1]
x
2
0
  d−1

   d+1
3 2
2
1 2
≥ 2Vd−1 [1]
2
d+1
2
  d−1



 d−1
3 2
2
1 2
≥ Vd−1 [1]
2
d+1
2


  d−1
2
3 2
.
≥ Vd−1 [1]
4
d+1

volume (Bd (x0 , 1) ∩ Bd (x1 , 1)) = 2

Using that

Vd−1

Vd−1 [1]
1
> √ , we divide the intersection volume through by Vd [1] to obtain,
Vd [1]
π
  d−1


3 2
2
√
4
π(d + 1)
  d2
1
3
≥
d+1 4

volume (Bd (x0 , 1) ∩ Bd (x1 , 1))
≥
volume (Bd (x0 , 1))

Lemma SM-G.7 (Packing balls in a ball). The number of non-intersecting balls of radius  which can be packed
d
into a ball of radius r in Rd is less than r
Proof. The technique used here is a loose version of that used in proving Lemma SM-G.5. The volume of Bd (0, )
d
is a factor (r/) smaller than that of Bd (0, r). As the balls of radius  are non-overlapping, the volume of their
union is simply the sum of their volumes. The result follow from the fact that the union of the balls of radius 
is contained within the ball of radius r.
Lemma SM-G.8 (Packing points in a ball). Given X ⊂ Bd (0, r) such that no two elements of X lie within a
d
distance of  of each other, |X | < 2r+
.

Proof. As no two elements lie within distance  of each other, balls of radius /2 centred at elements are nonintersecting. As each of the balls of radius /2 centred at elements of X lies entirely within Bd (0, r + /2), we
can apply Lemma (SM-G.7), arriving at the result.

SM-H

Pseudocode for trikmeds

In Alg. (6) we present trikmeds. It is decomposed into algorithms for initialisation (7), updating medoids (8),
assigning data to clusters (9) and updating bounds on the trimed derived bounds (10). Table 4 summarised all
of the variables used in trikmeds.
When there are no distance bounds, the location of the bottleneck in terms of distance calculations depends on
N/K 2 . If N/K  K, the bottleneck lies in updating medoids, which can be improved through the strategy
used in trimed. If N/K  K, the bottleneck lies in assigning elements to clusters, which is effectively handled
through the approach of Elkan (2003).

A Sub-Quadratic Exact Medoid Algorithm

Algorithm 6 trikmeds
initialise()
while not converged do
update-medoids()
assign-to-clusters()
update-sum-bounds()
end while

Algorithm 7 initialise
// Initialise medoid indices, uniform random sample without replacement (or otherwise)
{m(1), . . . , m(K)} ← uniform-no-replacement({1, . . . , N })
for k = 1 : K do
// Initialise medoid and set cluster count to zero
c(k) ← x(m(k))
v(k) ← 0
// Set sum of in-cluster distances to medoid to zero
s(k) ← 0
end for
for i = 1 : N do
for k = 1 : K do
// Tightly initialise lower bounds on data-to-medoid distances
lc (i, k) ← kx(i) − c(k)k
end for
// Set assignments and distances to nearest (assigned) medoid
a(i) ← arg mink∈{1,...,K} lc (i, k)
d(i) ← lc (i, a(i))
// Update cluster count
v(a(i)) ← v(a(i)) + 1
// Update sum of distances to medoid
s(a(i)) ← s(a(i)) + d(i)
// Initialise lower bound on sum of in-cluster distances to x(i) to zero
ls (i) ← 0
end for
V (0) ← 0
for k = 1 : K do
// Set cumulative cluster count
V (k) ← V (k − 1) + v(k)
// Initialise lower bound on in-cluster sum of distances to be tight for medoids
ls (m(k)) ← s(k)
end for
// Make clusters contiguous
contiguate()

James Newling, François Fleuret

Algorithm 8 update-medoids
for k = 1 : K do
for i = V (k − 1) : V (k) − 1 do
// If the bound test cannot exclude i as m(k)
if ls (i) < s(k) then
// Make ls (i) tight by computing and cumulating all in-cluster distances to x(i),
ls (i) ← 0
for i0 = V (k − 1) : V (k) − 1 do
˜ 0 ) ← kx(i) − x(i0 )k
d(i
˜ 0)
ls (i) ← ls (i) + d(i
end for
// Re-perform the test for i as candidate for m(k), now with exact sums. If i is the new best candidate,
update some cluster information
if ls (i) < s(k) then
s(k) ← ls (i)
m(k) ← i
for i0 = V (k − 1) : V (k) − 1 do
d(i0 ) ← kx(i) − x(i0 )k
end for
end if
// Use computed distances to i to improve lower bounds on sums for all samples in cluster k (see Figure
X)
for i0 = V (k − 1) : V (k) − 1 do
˜ 0 )v(k) − ls (i)|)
ls (i0 ) ← max (ls (i0 ), |d(i
end for
end if
end for
// If the medoid of cluster k has changed, update cluster information
if m(k) 6= V (k − 1) then
p(k) ← kc(k) − x(m(k))k
c(k) ← x(m(k))
end if
end for

A Sub-Quadratic Exact Medoid Algorithm

Algorithm 9 assign-to-clusters
// Reset variables monitoring cluster fluxes,
for k = 1 : K do
// the number of arrivals to cluster k,
∆n−in (k) ← 0
// the number of departures from cluster k,
∆n−out (k) ← 0
// the sum of distances to medoid k of samples which leave cluster k
∆s−out (k) ← 0
// the sum of distances to medoid k of samples which arrive in cluster k
∆s−in (k) ← 0
end for
for i = 1 : N do
// Update lower bounds on distances to medoids based on distances moved by medoids
for k = 1 : K do
l(i, k) = l(i, k) − p(k)
end for
// Use the exact distance of current assignment to keep bound tight (might save future calcs)
l(i, a(i)) = d(i)
// Record current assignment and distance
aold = a(i)
dold = d(i)
// Determine nearest medoid, using bounds to eliminate distance calculations
for k = 1 : K do
if l(i, k) < d(i) then
l(i, k) ← kx(i) − c(k)k
if l(i, k) < d(i) then
a(i) = k
d(i) = l(i, k)
end if
end if
end for
// If the assignment has changed, update statistics
if aold 6= a(i) then
v(aold ) = v(aold ) − 1
v(a(i)) = v(a(i)) + 1
ls (i) = 0
∆n−in (a(i)) = ∆n−in (a(i)) + 1
∆n−out (aold ) = ∆n−out (aold ) + 1
∆s−in (a(i)) = ∆s−in (a(i)) + d(i)
∆s−out (aold ) = ∆s−out (aold ) + dold
end if
end for
// Update cumulative cluster counts
for k = 1 : K do
V (k) ← V (k − 1) + v(k)
end for
contiguate()

James Newling, François Fleuret

Table 4: Table Of Notation For trikmeds
N
i
x(i)
K
k
m(k)
c(k)
n1 (i)
a(i)
d(i)
v(k)
V (k)
lc (i, k)
ls (i)
p(k)
s(k)

:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:

number of training samples
index of a sample, i ∈ {1, . . . , N }
sample i
number of clusters
index of a cluster, k ∈ {1, . . . , K}
index of current medoid of cluster k, m(k) ∈ {1, . . . , N }
current medoid of cluster k, that is c(k) = x(m(k))
cluster index of centroid nearest to x(i)
cluster to which x(i) is currently assigned
distance from x(i) to c(a(i))
number of samples assigned to cluster k
number of samples assigned to a cluster of index less than k + 1
lowerbound on P
distance from x(i) to m(k)
lowerbound on i0 :a(i0 )=a(i) kx(i0 ) − x(i)k
distance moved (teleported) by m(k) in last update
sum of distances of samples in cluster k to medoid k

Algorithm 10 update-sum-bounds
for k = 1 : K do
// Obtain absolute and net fluxes of energy and count, for cluster k
Jsabs (k) = ∆s−in (k) + ∆s−out (k)
Jsnet (k) = ∆s−in (k) − ∆s−out (k)
Jnabs (k) = ∆n−in (k) + ∆n−out (k)
Jnnet (k) = ∆n−in (k) − ∆n−out (k)
for i = V (k − 1) : V (k) − 1 do
// Update the lower bound on the sum of distances
ls (i) ← ls (i) − min(Jsabs (k) − Jnnet (k)d(i), Jnabs (k)d(i) − Jsnet (k))
end for
end for

SM-I

Datasets

• Birch1, Birch2 : Synthetic 2-D datasets available from https://cs.joensuu.fi/sipu/datasets/
• Europe : Border map of Europe available from https://cs.joensuu.fi/sipu/datasets/
• U-Sensor Net : Undirected 2-D graph data. Points drawn uniformly from
√ unit square, with an undirected
edge connecting points when the distance between them is less than 1.25 N
• D-Sensor Net : Directed 2-D graph data. Points drawn uniformly√from unit square, with directed edge
connecting points when the distance between them is less than 1.45 N , direction chosen at random.
• Europe rail : The European rail network, the shapefile is available at http://www.mapcruzin.com/
free-europe-arcgis-maps-shapefiles.htm. We extracted edges from the shapefile using networkx available at https://networkx.github.io/.
Algorithm 11 contiguate
// This function performs an in place rearrangement over of variables a, d, l, x and m
// The permutation applied to a, d, l and x has as result a sorting by cluster,
// a(i) = k if i ∈ {V (k − 1), V (k)} for k ∈ {1, . . . , K}
// and moreover that the first element of each cluster is the medoid,
// m(k) = V (k − 1) for k ∈ {1, . . . , K}

A Sub-Quadratic Exact Medoid Algorithm

• Pennsylvania road The road network of Pennsylvania, the edge list is available directly from https://snap.
stanford.edu/data/
• Gnutella Peer-to-peer network data, available from https://snap.stanford.edu/data/
• MNIST (0) The ‘0’s in the MNIST training dataset.
• Conflong The conflongdemo data is available from https://cs.joensuu.fi/sipu/datasets/
• Colormo The colormoments data is available at http://archive.ics.uci.edu/ml/datasets/Corel+
Image+Features
• MNIST50 The MNIST dataset, projected into 50-dimensions using a random projection matrix where each
of the 784 × 50 elements in the matrix is i.i.d. N (0, 1).
• S1, S2, S3, S4, A1, A2, A3 All of these synthetic datasets are available from https://cs.joensuu.fi/
sipu/datasets/.
• thyroid, yeast, wine, breast, spiral All of these real world datasets are available from https://cs.joensuu.
fi/sipu/datasets/.

SM-J

Scaling with dimension of TOPRANK and TOPRANK2

Recall the assumption (3) made for the TOPRANK and TOPRANK2 algorithms. The assumption states that as one
approaches the minimum energy E ∗ from above, the density of elements decreases. In other words, the lowest
energy elements stand out from the rest and are not bunched up with very similar energies.
Consider the case where elements are points in Rd . Suppose that the density fX of points around the medoid
is bounded by 0 < ρ0 ≤ fX ≤ ρ1 , and that the energy grows quadratically in radius about the medoid. Then,
as the number of points at radius  is O(d−1 ), the density (by energy) of points at radius  is O(d−2 ). Thus
for d = 1 the assumption for TOPRANK and TOPRANK does not hold, which results in poor performance for d = 1.
For d = 2, the assumption holds, as the density (by energy) of points is constant. For d ≥ 2, as d increases the
energy distribution becomes more and more favourable for TOPRANK and TOPRANK2, as the low ranking elements
become more and more distinct with low energies becoming less probable. This explains the observation that
TOPRANK scales well with dimension in Figure 3.

SM-K

Example where geometric median is a poor approximation of medoid

There is no guarantee that the geometric median is close to the set medoid. Moreover, the element in S
which is nearest to g(S) is not necessarily the medoid, as illustrated in the following example. Suppose S =
{x(1), . . . , x(20)} ⊂ R2 , with x(i) = (0, 1) for i ∈ {1, . . . , 9}, x(i) = (0, −1) for i ∈ {10, . . . , 18}, x(19) = (1/2, 0)
and x(20) = (−1/2, 0). The geometric
median is (0, 0) and the nearest points to the geometric
√
√ median, x(19)
and x(20) have energy 1 + 18 3/2 ≈ 16.6. However, points {x(1), . . . , x(18)} have energy 2 3/2 + 9 = 10.7.
Thus by choosing a point in S which is nearest to the geometric median, one is choosing the element with the
highest energy, the opposite of the medoid.
Note the above example appears to violate the assumptions required for O(N 3/2 ) convergence of trimed, as it
requires that the probability density function vanishes at the distribution median. Indeed, in Rd it is the case
that if the O(N 3/2 ) assumptions are satisfied, the set medoid converges to the geometric median, and so the
geomteric median is a good approximation. We stress however that the geometric median is only relevant in
vector spaces.

SM-L

Miscellaneous

Figure 7 illustrates the idea behind algorithm trimed, comments in the caption.

James Newling, François Fleuret

1N=1

N=4

N = 16

N = 64

N = 256

0
medoid
computed
excluded

-1
0

50

100

N

150

200

250

Figure 7: Eliminating samples as potential medoids using only type 1 elimination, where we assume that the
medoid and its energy E ∗ are known, and so the radius of the exclusion ball of an element x is E(x) − E ∗ .
Uniformly sampling from [−1, 1] × [−1, 1], energies are computed only if the sample drawn does not lie in the
exclusion zone (union of balls). If the energy at x is computed, the exclusion zone is augmented by adding
Bd (x, E(x) − E ∗ ). Top left to right: the distribution of samples which are computed and excluded. Bottom: the
1
times at which samples are computed. We prove that probability of computation at time n is O(n− 2 ).

