Energy Clustering
Guilherme França∗ and Joshua T. Vogelstein†
Johns Hopkins University

Abstract

arXiv:1710.09859v1 [stat.ML] 26 Oct 2017

Energy statistics was proposed by Székely in the 80’s inspired by the Newtonian gravitational
potential from classical mechanics, and it provides a hypothesis test for equality of distributions. It
was further generalized from Euclidean spaces to metric spaces of strong negative type, and more
recently, a connection with reproducing kernel Hilbert spaces (RKHS) was established. Here we
consider the clustering problem from an energy statistics theory perspective, providing a precise
mathematical formulation yielding a quadratically constrained quadratic program (QCQP) in the
associated RKHS, thus establishing the connection with kernel methods. We show that this QCQP
is equivalent to kernel k-means optimization problem once the kernel is fixed. These results imply
a first principles derivation of kernel k-means from energy statistics. However, energy statistics
fixes a family of standard kernels. Furthermore, we also consider a weighted version of energy
statistics, making connection to graph partitioning problems. To find local optimizers of such
QCQP we propose an iterative algorithm based on Hartigan’s method, which in this case has
the same computational cost as kernel k-means algorithm, based on Lloyd’s heuristic, but usually
with better clustering quality. We provide carefully designed numerical experiments showing the
superiority of the proposed method compared to kernel k-means, spectral clustering, standard
k-means, and Gaussian mixture models in a variety of settings.

∗

guifranca@gmail.com

†

jovo@jhu.edu

1

I.

INTRODUCTION

Energy statistics [1] is based on a notion of statistical potential energy between probability
distributions, in close analogy to Newton’s gravitational potential in classical mechanics. It
provides a model-free hypothesis test for equality of distributions which is achieved under
minimum energy. When probability distributions are different the statistical potential energy
diverges as sample size increases, while tends to a nondegenerate limit distribution when
probability distributions are equal. Energy statistics has been applied to several goodnessof-fit hypothesis tests, multi-sample tests of equality of distributions, analysis of variance
[2], nonlinear dependence tests through distance covariance and distance correlation, which
generalizes the Pearson correlation coefficient, and hierarchical clustering [3] by extending
Ward’s method of minimum variance. Moreover, in Euclidean spaces, an application of
energy statistics to clustering was already proposed [4]. We refer the reader to [1], and
references therein, for an overview of energy statistics theory and its applications.
In its original formulation, energy statistics has a compact representation in terms of
expectations of pairwise Euclidean distances, providing straightforward empirical estimates.
The notion of distance covariance was further generalized from Euclidean spaces to metric
spaces of strong negative type [5]. Furthermore, the missing link between energy distance
based tests and kernel based tests has been recently resolved [6], establishing an equivalence
between generalized energy distances to maximum mean discrepancies (MMD), which are
distances between embeddings of distributions in reproducing kernel Hilbert spaces (RKHS).
This equivalence immediately relates energy statistics to kernel methods often used in machine learning, and form the basis of our approach.
Clustering has such a long history in machine learning, making it impossible to mention
all important contributions in a short space. Perhaps, the most used method is k-means [7–
9], which is based on Lloyd’s heuristic [7] of assigning a data point to the cluster with closest
center. The only statistical information about each cluster comes from its mean, making it
sensitive to outliers. Nevertheless, k-means works very well when data is linearly separable
in Euclidean space. Gaussian mixture models (GMM) is another very common approach,
providing more flexibility than k-means, however, it still makes strong assumptions about
the distribution of the data.
To account for nonlinearities, kernel methods were introduced [10, 11]. A mercer kernel
2

[12] is used to implicitly map data points to a RKHS, then clustering can be performed in
the associated Hilbert space by using its inner product. However, the kernel choice remains
the biggest challenge since there is no principled theory to construct a kernel for a given
dataset, and usually a kernel introduces hyperparameters that need to be carefully chosen.
A well-known kernel based clustering method is kernel k-means, which is precisely k-means
formulated in the feature space [11]. Furthermore, kernel k-means algorithm [13, 14] is still
based on Loyd’s heuristic [7] of grouping points that are closer to a cluster center. We refer
the reader to [15] for a survey of clustering methods.
Although clustering from energy statistics, in Euclidean spaces, was considered in [4], the
precise optimization problem behind this approach remains elusive, as well as the connection with kernel methods. The main theoretical contribution of this paper is to fill this gap.
Since the statistical potential energy is minimum when distributions are equal, the principle
behind clustering is to maximize the statistical energy, enforcing probability distributions
associated to each cluster to be different from one another. We provide a precise mathematical formulation to this statement, leading to a quadratically constrained quadratic
program (QCQP) in the associated RKHS. This immediately establishes the connection between energy statistics based clustering, or energy clustering for short, with kernel methods.
Moreover, our formulation holds for general semimetric spaces of negative type. We also
show that such QCQP is equivalent to kernel k-means optimization problem, however, the
kernel is fixed by energy statistics. The equivalence between kernel k-means, spectral clustering, and graph partitioning problems is well-known [13, 14]. We further demonstrate how
these relations arise from a weighted version of energy statistics.
Our main algorithmic contribution is to use Hartigan’s method [16] to find local solutions
of the above mentioned QCQP, which is NP-hard in general. Hartigan’s method was also
used in [4], but without any connection to kernels. More importantly, the advantages of
Hartigan’s over Lloyd’s method was already demonstrated in some simple settings [17, 18],
but apparently this method did not receive the deserved attention. To the best of our
knowledge, Hartigan’s method was not previously employed together with kernel methods.
We provide a fully kernel based Hartigan’s algorithm for clustering, where the kernel is fixed
by energy statistics. We make clear the advantages of this proposal versus Lloyd’s method,
which kernel k-means is based upon and will also be used to solve our QCQP. We show that
both algorithms have the same time complexity, but Hartigan’s method in kernel spaces
3

offer several advantages. Furthermore, in the examples considered in this paper, it also
provides superior performance compared spectral clustering, which is more expensive and
in fact solves a relaxed version of our QCQP.
Our numerical results provide compelling evidence that Hartigan’s method applied to
energy clustering is more accurate and robust than kernel k-means algorithm. Furthermore,
our experiments illustrate the flexibility of energy clustering, showing that it is able to
perform accurately on data coming from very different distributions, contrary to k-means
and GMM for instance. More specifically, the proposed method performs closely to k-means
and GMM on normally distributed data, however, it is significantly better on data that is
not normally distributed. Its superiority in high dimensions is striking, being more accurate
than k-means and GMM even on Gaussian settings.

II.

BACKGROUND ON ENERGY STATISTICS AND RKHS

In this section we introduce the main concepts from energy statistics and its relation to
RKHS which form the basis of our work. For more details we refer the reader to [1] and [6].
iid

iid

Consider random variables in RD such that X, X 0 ∼ P and Y, Y 0 ∼ Q, where P and Q
are cumulative distribution functions with finite first moments. The quantity
E(P, Q) ≡ 2EkX − Y k − EkX − X 0 k − EkY − Y 0 k,

(1)

called energy distance [1], is rotationally invariant and nonnegative, E(P, Q) ≥ 0, where
equality to zero holds if and only if P = Q. Above, k · k denotes the Euclidean norm in RD .
Energy distance provides a characterization of equality of distributions, and E 1/2 is a metric

on the space of distributions.
The energy distance can be generalized as, for instance,
Eα (P, Q) ≡ 2EkX − Y kα − EkX − X 0 kα − EkY − Y 0 kα

(2)

where 0 < α ≤ 2. This quantity is also nonnegative, Eα (P, Q) ≥ 0. Furthermore, for
0 < α < 2 we have that Eα (P, Q) = 0 if and only if P = Q, while for α = 2 we have

E2 (P, Q) = 2kE(X) − E(Y )k2 which shows that equality to zero only requires equality of the
means, and thus E2 (P, Q) = 0 does not imply equality of distributions.
The energy distance can be even further generalized. Let X, Y ∈ X where X is an
arbitrary space endowed with a semimetric of negative type ρ : X × X → R, which is
4

required to satisfy
n
X

i,j=1

where Xi ∈ X and ci ∈ R such that

ci cj ρ(Xi , Xj ) ≤ 0,

Pn

i=1 ci

= 0. Then, X is called a space of negative type.

We can thus replace RD → X and kX − Y k → ρ(X, Y ) in the definition (1), obtaining the
generalized energy distance
E(P, Q) ≡ 2Eρ(X, Y ) − Eρ(X, X 0 ) − Eρ(Y, Y 0 ).

(3)

For spaces of negative type there exists a Hilbert space H and a map ϕ : X → H such that

ρ(X, Y ) = kϕ(X) − ϕ(Y )k2H . This allows us to compute quantities related to probability
distributions over X in the associated Hilbert space H. Even though the semimetric ρ may
not satisfy the triangle inequality, ρ1/2 does since it can be shown to be a proper metric. Our

energy clustering formulation, proposed in the next section, will be based on the generalized
energy distance (3).
There is an equivalence between energy distance, commonly used in statistics, and distances between embeddings of distributions in RKHS, commonly used in machine learning.
This equivalence was established in [6]. Let us first recall the definition of RKHS. Let H be
a Hilbert space of real-valued functions over X . A function K : X × X → R is a reproducing
kernel of H if it satisfies the following two conditions:
1. hx ≡ K(·, x) ∈ H for all x ∈ X .
2. hhx , f iH = f (x) for all x ∈ X and f ∈ H.
In other words, for any x ∈ X and any function f ∈ H, there is a unique hx ∈ H that
reproduces f (x) through the inner product of H. If such a kernel function K exists, then
H is called a RKHS. The above two properties immediately imply that K is symmetric
and positive definite. Indeed, notice that hhx , hy i = hy (x) = K(x, y), and by definition

hhx , hy i∗ = hhy , hx i, but since the inner product is real we have hhy , hx i = hhx , hy i, or
P
equivalently K(y, x) = K(x, y). Moreover, for any w ∈ H we can write w = ni=1 ci hxi
P
where {hxi }ni=1 is a basis of H. It follows that hw, wiH = ni,j=1 ci cj K(xi , xj ) ≥ 0, showing

that the kernel is positive definite. If G is a matrix with elements Gij = K(xi , xj ) this is
equivalent to G being positive semidefinite, i.e. v > G v ≥ 0 for any vector v ∈ Rn .
5

The Moore-Aronszajn theorem [19] establishes the converse of the above paragraph. For
every symmetric and positive definite function K : X ×X → R, there is an associated RKHS
HK with reproducing kernel K. The map ϕ : x 7→ hx ∈ HK is called the canonical feature
map. Given a kernel K, this theorem enables us to define an embedding of a probability
R
measure P into the RKHS as follows: P 7→ hP ∈ HK such that f (x)dP (x) = hf, hP i for
R
all f ∈ HK , or alternatively, hP ≡ K( · , x)dP (x). We can now introduce the notion of

distance between two probability measures using the inner product of HK , which is called

the maximum mean discrepancy (MMD) and is given by
γK (P, Q) ≡ khP − hQ kHK .

(4)

2
γK
(P, Q) = EK(X, X 0 ) + EK(Y, Y 0 ) − 2EK(X, Y )

(5)

This can also be written as [20]

iid

iid

where X, X 0 ∼ P and Y, Y 0 ∼ Q. From the equality between (4) and (5) we also have
hhP , hQ iHK = E K(X, Y ).
Thus, in practice, we can estimate the inner product between embedded distributions by
averaging the kernel function over sampled data.
The following important result shows that semimetrics of negative type and symmetric
positive definite kernels are closely related [21]. Let ρ : X × X → R and x0 ∈ X an arbitrary
but fixed point. Define
K(x, y) ≡ 21 [ρ(x, x0 ) + ρ(y, x0 ) − ρ(x, y)] .

(6)

Then, it can be shown that K is positive definite if and only if ρ is a semimetric of negative
type. We have a family of kernels, one for each choice of x0 . Conversely, if ρ is a semimetric
of negative type and K is a kernel in this family, then
ρ(x, y) = K(x, x) + K(y, y) − 2K(x, y)
= khx − hy k2HK

(7)

and the canonical feature map ϕ : x 7→ hx is injective [6]. When these conditions are satisfied
we say that the kernel K generates the semimetric ρ. If two different kernels generate the
same ρ they are said to be equivalent kernels.
6

Now we can state the equivalence between the generalized energy distance (3) and inner
products on RKHS, which is one of the main results of [6]. If ρ is a semimetric of negative
type and K a kernel that generates ρ, then replacing (7) into (3), and using (5), yields
2
(P, Q).
E(P, Q) = 2 [E K(X, X 0 ) + E K(Y, Y 0 ) − 2E K(X, Y )] = 2γK

Due to (4) we can compute the energy distance E(P, Q) between two probability distributions
using the inner product of HK .
Finally, let us recall the main formulas from generalized energy statistics for the test
statistic of equality of distributions [1]. Assume we have data X = {x1 , . . . , xn } where
S
xi ∈ X , and X is a space of negative type. Consider a disjoint partition X = kj=1 Cj ,

with Ci ∩ Cj = ∅. Each expectation in the generalized energy distance (3) can be computed
through the function
g(Ci , Cj ) ≡

1 XX
ρ(x, y),
ni nj x∈C y∈C
i

(8)

j

where ni = |Ci | is the number of elements in partition Ci . The within energy dispersion is
defined by
W ≡

k
X
nj
j=1

2

g(Cj , Cj ),

(9)

and the between-sample energy statistic is defined by
S≡
where n =

Pk

j=1

X

ni nj
[2g(Ci , Cj ) − g(Ci , Ci ) − g(Cj , Cj )] ,
2n
1≤i<j≤k

(10)

nj . Given a set of distributions {Pj }kj=1 , where x ∈ Cj if and only if x ∼ Pj ,

the quantity S provides a test statistic for equality of distributions [1]. When the sample size
is large enough, n → ∞, under the null hypothesis H0 : P1 = P2 = · · · = Pk we have that
S → 0, and under the alternative hypothesis H1 : Pi 6= Pj for at least two i 6= j, we have
that S → ∞. Note that this test does not make any assumptions about the distributions
Pj , thus it is said to be non-parametric or distribution-free.
One can make a physical analogy by thinking that points x ∈ Cj form a massive body
whose total mass is characterized by the distribution function Pj . The quantity S is thus a
potential energy of the from S(P1 , . . . , Pk ) which measures how different the distribution of
these masses are, and achieves the ground state S = 0 when all bodies have the same mass
distribution. The potential energy S increases as bodies have different mass distributions.
7

III.

CLUSTERING BASED ON ENERGY STATISTICS

This section contains the main theoretical results of this paper, where we formulate an
optimization problem for clustering based on energy statistics and RKHS introduced in the
previous section.
Due to the previous test statistic for equality of distributions, the obvious criterion for
clustering data is to maximize S which makes each cluster as different as possible from
the other ones. In other words, given a set of points coming from different probability
distributions, the test statistic S should attain a maximum when each point is correctly
classified as belonging to the cluster associated to its probability distribution. The following
straightforward result shows that maximizing S is, however, equivalent to minimizing W
which has a more convenient form.
Lemma 1. Let X = {x1 , . . . , xn } where each data point xi lives in a space X endowed with a
S
semimetric ρ : X × X → R of negative type. For a fixed integer k, the partition X = kj=1 Cj ,

where Ci ∩ Cj = ∅ for all i 6= j, maximizes the between-sample statistic S, defined in equation
(10), if and only if
min W (C1 , . . . , Ck ),

C1 ,...,Ck

(11)

where the within energy dispersion W is defined by (9).
Proof. From (9) and (10) we have

k
k 
k
X
1 X
1 X
nj ni g(Ci , Ci )
ni nj g(Ci , Cj ) +
n−
S+W =
2n i,j=1
2n i=1
j6=i=1
i6=j

=

k
1 X
1 XX
n
ni nj g(Ci , Cj ) =
ρ(x, y) = g(X, X).
2n i,j=1
2n x∈X y∈X
2

Note that the right hand side of this equation only depends on the pooled data, so it is a
constant independent of the choice of partition. Therefore, maximizing S over the choice of
partition is equivalent to minimizing W .
For a given k, the clustering problem amounts to finding the best partition of the data by
minimizing W . Notice that this is a hard clustering problem as partitions are disjoint. The
optimization problem (11) based on energy statistics was already proposed in [4]. However,
it is important to note that this is equivalent to maximizing S, which is the test statistic for
8

equality of distributions. In this current form, the relation with kernels and other clustering
methods is obscure. In the following, we show what is the explicit optimization problem
behind (11) in the corresponding RKHS, establishing the connection with kernel methods.
Based on the relation between kernels and semimetrics of negative type, assume that the
kernel K : X × X → R generates ρ. Define the Gram matrix


 K(x1 , x1 ) K(x1 , x2 ) · · · K(x1 , xn ) 


 K(x2 , x1 ) K(x2 , x2 ) · · · K(x2 , xn ) 
.
G≡
..
..
..


..
.


.
.
.


K(xn , x1 ) K(xn , x2 ) · · · K(xn , xn )

(12)

Let Z ∈ {0, 1}n×k be the label matrix, with only one nonvanishing entry per row, indicating

to which cluster (column) each point (row) belongs to. This matrix satisfies Z > Z = D,
where the diagonal matrix D = diag(n1 , . . . , nk ) contains the number of points in each

cluster. We also introduce the rescaled matrix Y ≡ ZD−1/2 . In component form they are
given by
Zij ≡



1 if xi ∈ Cj

Yij ≡


0 otherwise



 √1

nj


0

if xi ∈ Cj

.

(13)

otherwise

Throughout the paper, we use the notation Mi• to denote the ith row of a matrix M , and
M•j denotes its jth column. Our next result shows that the optimization problem (11) is
NP-hard since it is a quadratically constrained quadratic program (QCQP) in the RKHS.
Proposition 2. The optimization problem (11) is equivalent to
max Tr Y > G Y
Y



s.t. Y ≥ 0, Y > Y = I, Y Y > e = e,

(14)

where e = (1, 1, . . . , 1)> ∈ Rn is the all-ones vector, and G is the Gram matrix (12).
Proof. From (7), (8), and (9) we have

k
k X
X
1X 1 X
1 X
W (C1 , . . . , Ck ) =
K(x, x) −
ρ(x, y) =
K(x, y) .
2 j=1 nj x,y∈C
nj y∈C
j=1 x∈C
j

j

(15)

j

Note that the first term is global so it does not contribute to the optimization problem.
Therefore, minimizing (15) is equivalent to
k
X
1 X
max
K(x, y).
C1 ,...,Ck
n
j=1 j x,y∈C
j

9

(16)

But

X

K(x, y) =

n X
n
X

Zpj Zqj Gpq = (Z > G Z)jj ,

p=1 q=1

x,y∈Cj

−1
where we used the definitions (12) and (13). Notice that n−1
= Djj
, where the diagonal
j

matrix D = diag(n1 , . . . , nk ) contains the number of points in each cluster, thus the objective


P
−1
function in (16) is equal to kj=1 Djj
Z > GZ jj = Tr D−1 Z > GZ . Now we can use the

cyclic property of the trace, and by the definition of the matrix Z in (13), we obtain the
following integer programing problem:
max Tr
Z



ZD−1/2

>

G ZD−1/2



s.t. Zij ∈ {0, 1},

Pk

j=1

Zij = 1,

Pn

i=1

Zij = nj .

(17)

Now we write this in terms of the matrix Y = ZD−1/2 . The objective function immedi
ately becomes Tr Y > G Y . Notice that the above constraints imply that Z T Z = D, which

in turn gives D−1/2 Y T Y D−1/2 = D, or Y > Y = I. Also, every entry of Y is positive by
definition, Y ≥ 0. Now it only remains to show the last constraint in (14), which comes

from the last constraint in (17). In matrix form this reads Z T e = De. Replacing Z = Y D1/2

we have Y > e = D1/2 e. Multiplying this last equation on the left by Y , and noticing that
Y D1/2 e = Ze = e, we finally obtain Y Y > e = e. Therefore, the optimization problem (17)
is equivalent to (14) .
Based on Proposition 2, to group data X = {x1 , . . . , xn } into k clusters we first compute

the Gram matrix G and then solve the optimization problem (14) for Y ∈ Rn×k . The
ith row of Y will contain a single nonzero element in some jth column, indicating that
xi ∈ Cj . This optimization problem is nonconvex, and also NP-hard, thus a direct approach
is computational prohibitive even for small datasets. However, one can find approximate
solutions by relaxing some of the constraints, or obtaining a relaxed SDP version of it. For
instance, the relaxed problem
max Tr Y > G Y
Y



s.t. Y > Y = I

has a well-known closed form solution Y ? = U R, where the columns of U ∈ Rn×k contain
the top k eigenvectors of G corresponding to the k largest eigenvalues λ1 ≥ λ2 ≥ . . . ≥ λk ,

and R ∈ Rk×k is an arbitrary orthogonal matrix. The resulting optimal objective function
 P
assumes the value max Tr Y ? > G Y ? = ki=1 λi . Spectral clustering is based on the above
10

approach, where one further normalize the rows of Y ? , then cluster the resulting rows as
data points. A procedure on these lines was proposed in the seminal papers [22, 23].
Note that the optimization problem (14) based on energy statistics is valid for data living
in an arbitrary space of negative type, where a semimetric ρ, and thus the kernel K, are
assumed to be known. The standard energy distance (2) fixes a family of choices in Euclidean
spaces given by
ρα (x, y) = kx − ykα ,

Kα (x, y) = 21 (kxkα + kykα − kx − ykα ) ,

for 0 < α ≤ 2 and we fix x0 = 0 in (6). The same would be valid for data living in a more
general semimetric space (X , ρ) where ρ fixes the kernel. In practice, the clustering quality
strongly depend on the choice of a suitable ρ. Nevertheless, if prior information is available
to make this choice, it can be immediately incorporated in the optimization problem (14).

Relation to Kernel k-Means

One may wonder how energy clustering relates to the well-known kernel k-means problem1
which is extensively used in machine learning. For a positive semidefinite Gram matrix G,
as defined in (12), there exists a map ϕ : X → HK such that K(x, y) = hϕ(x), ϕ(y)i. Kernel
k-means optimization problem, in feature space, is defined by


k X
X
2
min J(C1 , . . . , Ck ) ≡
kϕ(x) − ϕ(µj )k
C1 ,...,Ck

where µj =

1
nj

P

x∈Cj

(18)

j=1 x∈Cj

x is the mean of cluster Cj in the ambient space. Notice that the above

objective function is strongly tied to the idea of minimizing distances between points and
cluster centers, which arises from k-means objective function based on Lloyd’s method [7].
It is known [13, 14] that problem (18) can be cast into a trace maximization in the same
form as (14). The next result makes this explicit, showing that (11) and (18) are actually
equivalent.
Proposition 3. For a fixed kernel, the clustering optimization problem (11) based on energy statistics is equivalent to the kernel k-means optimization problem (18), and both are
equivalent to (14).
1

When we refer to kernel k-means problem we mean specifically the optimization problem (18), which
should not be confused with kernel k-means algorithm that is just one possible recipe to solve (18). The
distinction should also be clear from the context.

11

Proof. Notice that kϕ(x)−ϕ(µj )k2 = hϕ(x), ϕ(x)i−2hϕ(x), ϕ(µj )i+hϕ(µj ), ϕ(µj )i, therefore

k X
X
1 X
2 X
K(x, y) + 2
K(y, z) .
J=
K(x, x) −
nj y∈C
nj y,z∈C
j=1 x∈C
j

j

(19)

j

The first term is global so it does not contribute to the optimization problem. Notice that
P
P
P
the third term gives x∈Cj n12 y,z∈Cj K(y, z) = n1j y,z∈Cj K(y, z), which is the same as the
j

second term. Thus, problem (18) is equivalent to

k
X
1 X
K(x, y)
C1 ,...,Ck
n
j
j=1
x,y∈C

max

j

which is exactly the same as (16) from the energy statistics formulation. Therefore, once the
kernel K is fixed, the function W given by (9) is the same as J in (18). The remaining of the
proof proceeds as already shown in the proof of Proposition 2, leading to the optimization
problem (14).
The above result shows that kernel k-means optimization problem is equivalent to the
clustering problem formulated in the energy statistics framework, when operating on the
same kernel. Notice, however, that energy statistics is valid for arbitrary semimetric spaces
of negative type, fixing the kernel function in the associated RKHS, which is guaranteed to
be positive definite. On the other hand, kernel k-means (18) by itself is just an heuristic
approach that does not make any explicit mention to the kernel. Based on Proposition 3
one may view kernel k-means as being derived from the energy statistics framework.
Kernel k-means, spectral clustering, and graph partitioning problems such as ratio association, ratio cut, and normalized cut are all equivalent to a QCQP of the form (14)
[13, 14]. One can thus use kernel k-means algorithm to solve these problems as well. This
correspondence involves a weighted version of problem (14), that will be demonstrated in
the following from the perspective of energy statistics.

IV.

CLUSTERING BASED ON WEIGHTED ENERGY STATISTICS

We now generalize energy statistics to incorporate weights associated to each data point.
Let w(x) be a weight function associated to point x ∈ X . Define
g(Ci , Cj ) ≡

1 XX
w(x)w(y)ρ(x, y),
si sj x∈C y∈C
i

j

12

si ≡

X

x∈Ci

w(x).

(20)

Replace this function in the formulas (9) and (10), with ni → si and n → s, where s =
Pk
j=1 sj . With these changes Proposition 1 remains the unaltered, so the clustering problem
becomes



k
X
sj
g(Cj , Cj )
min W (C1 , . . . , Ck ) ≡
C1 ,...,Ck
2
j=1

where now g is given by (20). Define the following matrices and vector:


 √1
if xi ∈ Cj
sj
Yij ≡
,
W ≡ diag(w1 , . . . , wn ),
H ≡ W 1/2 Y,

0
otherwise

(21)

ω ≡ We, (22)

where wi = w(xi ) and e ∈ Rn is the all-ones vector.The analogous of Proposition 2 is as
follows.
Proposition 4. The weighted energy clustering given by problem (21) is equivalent to

max Tr H > (W 1/2 GW 1/2 )H

s.t. H ≥ 0, H > H = I, HH > ω = ω,

H

(23)

where G is the Gram matrix (12), ω = (w1 , . . . , wn )T contains the weights of each point, and
W = diag(ω).
Proof. Replacing (7) and eliminating the global terms which do not contribute, the optimization problem (21) becomes
k
X
1 XX
w(x)w(y)K(x, y).
C1 ,...,Ck
s
j
j=1
x∈C y∈C

max

j

j

This objective function can be written as
√
k
n
n
k X
n X
n
>√
X
X
Zjp
wp 1/2
wq Zqj
1 XX
1/2
wp wq Zpj Zqj Gpq =
wp Gpq wq
√
√
s
sj
sj
j=1 j p=1 q=1
j=1 p=1 q=1
=

k
X
j=1

H > W 1/2 GW 1/2 H



jj


= Tr H > W 1/2 GW 1/2 H .

To obtain the constraints, note that Hij ≥ 0 by definition, and
(H > H)ij =

n
n
X
1
δij X
w` Z`i Z`j =
w` Z`i = δij ,
Y`i W`` Y`j = √ √
s
s
i sj
i
`=1
`=1
`=1

n
X

13

where δij = 1 if i = j and δij = 0 if i 6= j is the Kronecker delta. Therefore, H > H = I. This
is a constraint on the rows of H. To obtain a condition on its columns observe that
√
wp wq

k

if both xp , xq ∈ Ci
X

Z
Z
√
pj qj
si
=
H > H pq = wp wq

sj
0
j=1
otherwise.

Therefore, (H > HW 1/2 )pq =

√

wp wq s−1
if both points xp and xq belong to the same cluster,
i

which we denote by Ci for some i ∈ {1, . . . , k}, and (H > HW 1/2 )pq = 0 otherwise. Thus, the
pth line of this matrix is nonzero only on entries corresponding to points that are in the same
Pn
√
√
wp ,
cluster as xp . If we sum over the columns of this line we obtain wp s−1
i
q=1 wq Zqi =
or equivalently HH > W 1/2 e = W 1/2 e, which gives the constraint HH > ω = ω.
Connection with Graph Partitioning

The relation between kernel k-means and graph partitioning problems is known [13, 14].
For conciseness, we repeat a similar analysis due to the relation of these problems to energy
statistics and RKHS, which provides a different perspective.
Consider a graph G = (V, E, A) where V is the set of vertices, E the set of edges, and A
is an affinity matrix of the graph, which measures the similarities between pairs of nodes.
Thus, Aij 6= 0 if (i, j) ∈ E, and Aij = 0 otherwise. We also associate weights to every vertex,
P
wi = w(i) for i ∈ V, and let sj = i∈Cj wi , where Cj ⊆ V is one partition of V. Let
links(C` , Cm ) ≡

X

Aij .

i∈C` ,j∈Cm

We want to partition the set of vertices V into k disjoint subsets, V =
generalized ratio association problem is given by
max

Ci ,...,Ck

Sk

j=1

k
X
links(Cj , Cj )
j=1

sj

Cj . The

(24)

and maximizes the within cluster association. The generalized ratio cut problem
min

Ci ,...,Ck

k
X
links(Cj , V\Cj )

sj

j=1

(25)

minimizes the cut between clusters. These two problems are equivalent, in analogous way
as minimizing (9) is equivalent to maximizing (10) as shown in Proposition 1. Here this is
14

due to the equality links(Cj , V\Cj ) = links(Cj , V) − links(Cj , Cj ). Several graph partitioning
methods [22, 24–26] can be seen as a particular case of (24) or (25).
Consider the ratio association problem (24), whose objective function can be written as
k X
k
n X
n
>
X
X

Zjp
Zqj
1 XX
>
Apq =
√ Apq √ = Tr Y AY ,
s
sj
sj
j=1 p=1 q=1
j=1 j p∈C q∈C
j

j

with Z defined in (13) and Y in (22). Therefore, the ratio association problem can be written
in the form (23), i.e.
max Tr H > W −1/2 AW −1/2 H
H



s.t. H ≥ 0, H > H = I, HH > ω = ω.

This is exactly the same problem as weighted energy clustering with G = W −1 AW −1 .
Assuming this matrix is positive semidefinite, this generates a semimetric (7) for graphs
given by
ρ(i, j) =

Aii Ajj
2Aij
+ 2 −
2
wi
wj
wi wj

or

ρ(i, j) = −

2Aij
wi wj

(26)

for vertices i, j ∈ V, and where in the second equation we assume the graph has no self-loops,
i.e. Aii = 0. Using (26) in the energy statistics formulation allows one to make inference on
graphs. Above, the weight wi = w(i) of node i ∈ V can be, for instance, its degree wi = d(i).
V.

TWO-CLASS PROBLEM IN ONE DIMENSION

Before stating a general algorithm to solve the optimization problem (14) we first consider
the simplest possible case which is one-dimensional data and a two-class problem. This will
be useful to test energy clustering on a simple setting.
Fixing ρ(x, y) = |x − y| according to the standard energy distance, we can actually
compute the function (8) in O(n log n) and minimize W directly. This is done by noting
that

|x − y| = (x − y)1x≥y − (x − y)1x<y

= x (1x≥y − 1x<y ) + y (1y>x − 1y≤x )

where we have the indicator function defined by 1A = 1 if A is true, and 1A = 0 otherwise.
Let C be a partition with n elements. Using the above distance we have
g (C, C) =

1 XX
x (1x≥y + 1y>x − 1x≥y − 1x<y ) .
n2 x∈C y∈C
15

Algorithm 1 E 1D -clustering algorithm to find local solutions to the optimization problem
(11) for a two-class problem in one dimension.
input data X
output label matrix Z
e = [x1 , . . . , xn ]
1: sort X obtaining X
2: for j ∈ [1, . . . , n] do
3:
4:

Ce1,j ← [xi : i = 1, . . . , j], and Ce2,j ← [xi : i = j + 1, . . . , n]

W (j) ← W Ce1,j , Ce2,j , from (27)

5: end for

6: j ? ← arg minj W (j)
7: Zj• ← (1, 0) if j ≤ j ? , and Zj• ← (0, 1) otherwise, for j = 1, . . . , n

The sum over y can be eliminated since each term in the parenthesis is simply counting the
number of elements in C that satisfy the condition of the indicator function. Assuming that
we first order the data in C, obtaining Ce = [xj ∈ C : x1 ≤ x2 ≤ · · · ≤ xn ], we get
n

2 X
e
e
(2` − 1 − n)x` .
g C, C = 2
n `=1


e Ce is O(n) and the cost of sorting the data is at the
Note that the cost of computing g C,
S
most O(n log n). Assuming that each partition is ordered, X = kj=1 Cej , the within energy
dispersion can be written explicitly as



W Ce1 , . . . , Cek =

nj
k X
X
2` − 1 − nj
j=1 `=1

nj

x` .

(27)

For a two-class problem we can use the formula (27) to cluster the data through a simple
e Then we compute (27)
algorithm as follows. We first order the entire dataset, X → X.

e and pick the point which gives the minimum value of W . This
for each possible split of X

procedure is described in Algorithm 1 and called E 1D -clustering. Note that this algorithm
is deterministic, however, it only works for one-dimensional data with Euclidean distance.
The total complexity of E 1D -clustering is O(n log n + n2 ) = O(n2 ).
Assuming the true label matrix Z is available, a direct measure of how different the
16

estimated matrix Ẑ is from Z, up to label permutations, is given by
n

k

1 XX
accuracy(Ẑ) ≡ max
Ẑiσ(j) Zij
σ n
i=1 j=1

(28)

where σ is a permutation of the k cluster groups. The accuracy is always between [0, 1],
where 1 corresponds to all points correctly clustered, and 0 to all points wrongly clustered.
For a balanced two-class problem the value 1/2 correspond to chance.
We now consider two simple experiments where we sample n points from a two-class
mixture. We plot the average accuracy (28) versus n, with error bars indicating standard
error. The data is clustered using E 1D -clustering algorithm, GMM and k-means. For GMM
and k-means we use the implementations from the well-known scikit-learn library in Python
[27], where k-means is initialized through k-means++ procedure [28], and GMM is initialized
with the output of k-means. We run both algorithms 5 times with different initializations
and pick the answer with best objective function value. Notice that E 1D -clustering does not
require random initialization so we only run it once. For each n we use use 100 Monte Carlo
runs. In Fig. 1a we have the results for data sampled from the Gaussian mixture


iid
µ1 = 1.5, σ1 = 0.3, µ2 = 0, σ2 = 1.5.
x ∼ 21 N µ1 , σ12 + 12 N µ2 , σ22 ,

(29)

In this case the optimal accuracy obtained from Bayes classification error is ≈ 0.88, indicated
by the dashed line in the plot. The three methods perform closely, with a slight advantage
of GMM, as expected since it is a consistent model to the data, and E 1D -clustering performs
slightly better than k-means. In Fig. 1c we show a density estimation from clustering
1000 points from this mixture using the three algorithms. Notice that all of them are able
to distinguish the two classes. On the other hand, in Fig. 1b we consider a mixture of
lognormal distributions,




iid
x ∼ 21 exp N µ1 , σ12 + 12 exp N µ2 , σ22 ,

µ1 = 1.5, σ1 = 0.3, µ2 = 0, σ2 = 1.5. (30)

The optimal Bayes accuracy is again ≈ 0.88. We can now see that E 1D -clustering is still
very accurate, while GMM and k-means basically cluster at chance. Density estimation
after clustering 1000 points this mixture using the three algorithms are is shown in Fig. 1d.
Note that only E 1D -clustering was able to distinguish the two classes. k-means and GMM
put most of the points in a single cluster, and points on the tail of the second component of
(30) in the other cluster. The experiments of Fig. 1 illustrate how energy clustering is more
flexible compared to k-means and GMM.
17

0.90

0.90

0.80

0.85

accuracy

accuracy

0.85
1D

E -clustering
k-means
GMM

0.80

0.75

E 1D -clustering
k-means
GMM

0.70
0.65
0.60
0.55

0.75
100

200

300

400

500

600

700

800

0.50

100

200

300

0.8

truth
E 1D
k-means
GMM

700

800

0.5
0.4
0.3

0.4

0.2

0.2

0.1
−4

truth
E 1D
k-means
GMM

0.6

0.6

0.0
−6

600

0.7

1.4

1.0

500

(b)

(a)
1.2

400

# points

# points

−2

0

2

4

6

x

0.0
−2

0

2

4

6

8

10

12

14

x

(c)

(d)

FIG. 1. E 1D -clustering versus k-means and GMM. (a,b) We plot the mean accuracy (28) over
100 Monte Carlo trials, versus the number of sampled points. Error bars are standard error. The
dashed line indicates Bayes accuracy (≈ 0.88 in both cases). (a) Clustering results for data normally
distributed as in (29). (b) Data lognormally distributed as in (30). (c) Density estimation of each
component in the mixture (29) after clustering 1000 sampled points using the three algorithms,
compared to the ground truth. (d) The same but for lognormal data (30).
VI.

ITERATIVE ALGORITHMS FOR ENERGY CLUSTERING

In this section we introduce an iterative algorithm to find a local maximizer of the optimization problem (14). Due to Proposition 3 we can also find an approximate solution
by the well-known kernel k-means algorithm based on Lloyd’s heuristic [13, 14], which for
convenience will also be restated in the present context.
Consider the optimization problem (16) written as
max

{C1 ,...,Ck }




k
X
Qj
Q=
,
nj
j=1
18

Qj ≡

X

x,y∈Cj

K(x, y),

(31)

where Qj represents an internal energy cost of cluster Cj , and Q is the total energy cost
where each Qj is weighted by the inverse of the number of points in Cj . For a data point xi
we denote its own energy cost with the entire cluster C` by
Q` (xi ) ≡

X

y∈C`

K(xi , y) = Gi• · Z•` ,

where we recall that Gi• (G•i ) denotes the ith row (column) of matrix G.

Lloyd’s Method for Energy Clustering

To optimize kernel k-means objective function (19) we remove the global term and define
the function
J (`) (xi ) ≡

2
1
Q − Q` (xi ).
2 `
n`
n`

(32)

We are thus solving
min
Z

n X
k
X

Zi` J (`) (xi ).

i=1 `=1

One possible strategy is to assign xi to cluster Cj ? according to
j ? = arg min J (`) (xi ).
`=1,...,k

This is done for every data point xi and repeated until convergence, i.e. until no new
assignments are made. The entire procedure is described in Algorithm 2, which we name
E L -clustering to emphasize that we are optimizing the within energy function W based on
Lloyd’s method [7]. It can be shown that this algorithm converges provided G is positive
semidefinite.
E L -clustering is precisely kernel k-means algorithm [13, 14] but written more concisely and
with the kernel induced by energy statistics. Indeed, recalling that K(x, y) = hϕ(x), ϕ(y)i
where ϕ : X → HK is the feature map, we have from (32) that
J (`) (xi ) = hϕ(µ` ), ϕ(µ` )i − 2hϕ(xi ), ϕ(µ` )i = kϕ(xi ) − ϕ(µ` )k2 − kϕ(xi )k2 ,
where µ` =

1
n`

P

x∈C`

x is the mean of cluster C` . Therefore, min` J (`) (xi ) = min` kϕ(xi ) −

ϕ(µ` )k2 , i.e. we are assigning xi to the cluster with closest center (in feature space), which
is the familiar Lloyd’s heuristic approach that kernel k-means is based upon.
19

Algorithm 2 E L -clustering is Lloyd’s method for energy clustering, which is precisely
kernel k-means algorithm, with the kernel induced by energy statistics. This procedure
finds local solutions to the optimization problem (14).
input number of clusters k, Gram matrix G, initial label matrix Z ← Z0
output label matrix Z
1: q ← (Q1 , . . . , Qk )> have the costs of each cluster, defined in (31)
2: n ← (n1 , . . . , nk )> have the number of points in each cluster
3: repeat
4:

for i = 1, . . . , n do

5:

let j be such that xi ∈ Cj

6:

j ? ← arg min`=1,...,k J (`) (xi ), where J (`) (xi ) is defined in (32)

7:

if j ? 6= j then

8:

move xi to Cj ? : Zij ← 0 and Zij ? ← 1

9:

update n: nj ← nj − 1 and nj ? ← nj ? + 1

10:

update q: qj ← qj − 2Qj (xi ) and qj ? ← qj ? + 2Qj ? (xi )

11:

end if

12:

end for

13: until convergence

To check the complexity of E L -clustering, notice that to compute the second term of

J (`) (xi ) in (32) requires O(n` ) operations, and although the first term requires O(n2` ) it

only needs to be computed once outside loop through data points (step 1 of Algorithm 2).
Therefore, the time complexity of E L -clustering is O(nk max` n` ) = O(kn2 ). For a sparse
Gram matrix G having n0 nonzero elements this complexity can be further reduced to O(kn0 ).

Hartigan’s Method for Energy Clustering

We now consider Hartigan’s method [16] applied to the optimization problem in the
form (31), which gives a local solution to the QCQP defined in (14). The method is based
in computing the maximum change in the total cost function Q when moving each data
20

point to another cluster. More specifically, suppose point xi is currently assigned to cluster
Cj yielding a total cost function denoted by Q(j) . Moving xi to cluster C` yields another

total cost function denoted by Q(`) . We are interested in computing the maximum change
∆Qj→` (xi ) ≡ Q(`) − Q(j) , for ` 6= j. From (31), by explicitly writing the costs related to
these two cluster we obtain
−

∆Q

j→`

Qj
Qj
Q`
Q+
`
+
−
−
(xi ) =
n` + 1 nj − 1
nj
n`

−
where Q+
` denote the cost of the new `th cluster with the point xi added to it, and Qj is

the cost of new jth cluster with xi removed from it. Noting that Q+
` = Q` + 2Q` (xi ) + Gii
and Q−
j = Qj − 2Qj (xi ) + Gii , we get the formula




Qj
1
Q`
1
j→`
− 2Qj (xi ) + Gii −
− 2Q` (xi ) − Gii .
∆Q (xi ) =
nj − 1 nj
n` + 1 n`

(33)

Therefore, if ∆Qj→` (xi ) > 0 we get closer to a maximum of (31) by moving xi to C` , otherwise
we keep xi in Cj .
We thus propose the following algorithm. We start with an initial configuration for the
label matrix Z, then for each point xi we compute the cost of moving it to another cluster
C` , i.e. ∆Qj→` (xi ) for ` = 1, . . . , k with ` 6= j, where j denotes the index of its current
partition, x ∈ Cj . Hence, we choose
j ? = arg max ∆j→` (xi ).
`=1,...,k | `6=j
?

If ∆Qj→j (xi ) > 0 we move xi to cluster Cj ? , otherwise we keep xi in its original cluster Cj .
This process is repeated until no points are assigned to new clusters. The entire procedure
is explicitly described in Algorithm 3, which we denote E H -clustering to emphasize that it is
based on Hartigan’s method. This method automatically ensures that the objective function
is monotonically increasing at each iteration, and consequently the algorithm converges in
a finite number of steps.
The complexity analysis of E H -clustering is the following. Computing the Gram matrix

G requires O(Dn2 ) operations, where D is the dimension of each data point and n is the

data size. However, both algorithms E L - and E H -clustering assume that G is given. There
are more efficient methods to compute G, specially if it is sparse, but we will not consider
this further and just assume that G is given. The computation of each cluster cost Qj has
complexity O(n2j ), and overall to compute q we have O(n21 + · · · + n2k ) = O(k maxj n2j ).
21

Algorithm 3 E H -clustering is Hartigan’s method for energy clustering. This algorithm
finds local solutions to the optimization problem (14). The steps 6 and 10 are different
than E L -clustering described in Algorithm 2.

input number of clusters k, Gram matrix G, initial label matrix Z ← Z0
output label matrix Z
1: q ← (Q1 , . . . , Qk )> have the energy costs of each cluster, defined in (31)
2: n ← (n1 , . . . , nk )> have the number of points in each cluster
3: repeat
4:

for i = 1, . . . , n do

5:

let j be such that xi ∈ Cj

6:

j ? ← arg max`=1,...,k | `6=j ∆Qj→` (xi ) using (33)

7:

if ∆Qj→j (xi ) > 0 then

?

8:

move xi to Cj ? : Zij ← 0 and Zij ? ← 1

9:

update n: nj ← nj − 1 and nj ? ← nj ? + 1

10:

update q: qj ← qj − 2Qj (xi ) + Gii and qj ? ← qj ? + 2Qj ? (xi ) + Gii

11:

end if

12:

end for

13: until convergence

These operations only need to be performed a single time. For each point xi we need to
compute Qj (xi ) once, which is O(nj ), and we need to compute Q` (xi ) for each ` 6= j. The
cost of computing Q` (xi ) is O(n` ), thus the cost of step 6 in Algorithm 3 is O(k max` n` ) for
` = 1, . . . , k. For the entire dataset this gives a time complexity of O(nk max` n` ) = O(kn2 ).

Note that this is the same cost as in E L -clustering, or kernel k-means algorithm. Again, if
G is sparse this can be reduced to O(kn0 ) where n0 is the number of nonzero entries of G.
In the following we mention some important known results about Hartigan’s method.
Theorem 5 (Telgarsky-Vattani [17]). Hartigan’s method has the cost function strictly decreasing in each iteration. Moreover, if n > k then
1. the resulting partition has no empty clusters, and
22

2. the resulting partition has distinct means.
Neither of these two conditions are guaranteed to be satisfied by Lloyd’s method, and
consequently by E L -clustering algorithm. The next result indicates that Hartigan’s method
can potentially escape local optima of Lloyd’s method.
Theorem 6 (Telgarsky-Vattani [17]). The set of local optima of Hartigan’s method is a
(possibly strict) subset of local optima of Lloyd’s method.
The above theorem implies that E L -clustering cannot improve on a local optima of E H -

clustering. On the other hand, E H might improve on a local optima of E L . Lloyd’s method

forms Voronoi partitions, while Hartigan’s method groups data in regions formed by the
intersection of spheres called circlonoi cells. It can be shown that the circlonoi cells are
contained within a smaller volume of a Voronoi cell, and this excess volume grows exponentially with the dimension of X [17, Theorems 2.4 and 3.1]. Points in this excess volume
force Hartigan’s method to iterate, contrary to Lloyd’s method. Therefore, Hartigan’s can
escape local optima of Lloyd’s. Moreover, this improvement should be more prominent as
dimension increases. Also, the improvement grows as the number of clusters k increases.
The empirical results of [17] show that an implementation of Hartigan’s method has comparable execution time to an implementation of Lloyd’s method, but no explicit complexity
was provided. We show that both E L - and E H -clustering have the same time complexity. To
the best of our knowledge, Hartigan’s method was not previously considered together with
kernels, as we are proposing in E H -clustering algorithm.
In [18], Hartigan’s method was applied to k-means problem with any Bregman divergence.
It was shown that the number of Hartigan’s local optima is upper bounded by O(1/k) [18,
Proposition 5.1]. In addition, it was provided examples where any initial partition correspond to a local optima of Lloyd’s method, while the number of local optima in Hartigan’s
method is small and correspond to true partitions of the data. Empirically, the number of
Hartigan’s local optima was considerably smaller than the number of Lloyd’s local optima.
The above results indicate that Hartigan’s method provides several advantages over
Lloyd’s method, a fact that will also be supported by our numerical experiments in the
next section where E H outperforms of E L (kernel k-means) in several settings, specially in
high dimensions.
23

VII.

NUMERICAL EXPERIMENTS

The main goal of this section is threefold. First, we want to compare E H -clustering in

Euclidean space to k-means and GMM. Second, we want to compare E H -clustering, based

on Hartigan’s method, to E L -clustering or kernel k-means, based on Lloyd’s method, and
also to spectral clustering, when they all operate on the same kernel. Third, we want to
illustrate the flexibility provided by energy clustering, which is able to cluster accurately in
different settings while keeping the same kernel.
The following experimental setup holds unless specified otherwise. We consider E H -

clustering, E L -clustering and spectral clustering with the following semimetrics and corresponding generating kernels:
ρα (x, y) = kx − ykα ,

Kα (x, y) = 21 (kxkα + kykα − kx − ykα ) ,

(34)

ρeσ (x, y) = 2 − 2e−

e σ (x, y) = e−
K

kx−yk
2σ ,

(35)

kx−yk2
2σ 2 .

(36)

ρbσ (x, y) = 2 − 2e−

kx−yk
2σ ,
kx−yk2
2σ 2 ,

b σ (x, y) = e−
K

The relation between kernel and semimetric is given by formula (6) where we fix x0 = 0. The
standard ρ1 , from the original energy distance (1), will always be present in the experiments
as a reference, being the implied choice unless explicitly mentioned. For k-means, GMM
and spectral clustering we use the robust implementations of scikit-learn library [27], where
k-means is initialized with k-means++ [28], and GMM with the output of k-means, making
it more robust and preventing it from breaking in high dimensions. Spectral clustering
implementation is based on [22]. We implemented E L -clustering as described in Algorithm 2,

and E H -clustering as described in Algorithm 3. Both will also be initialized with k-means++.
We run the algorithms 5 times with different initializations, picking the result with best
objective function value. We evaluate clustering quality by the accuracy (28) based on the
true labels. For each setting we show the average accuracy over 100 Monte Carlo trials, with
error bars indicating standard error.
We briefly mention that we compared E H -clustering, as described in Algorithm 3, to

E 1D -clustering, described in Algorithm 1, for several univariate distributions. Both perform
very closely. However, we omit these results since we will analyse more interesting scenarios
in high dimensions.
From the results of [17], summarized in the end of the previous section, we expect the
24

x1

−3

5
4
3
2
1
0
−1
−2
−3

x2

4
3
2
1
0
−1
−2
−3
−4

x2

8
6
4
2
0
−2
−4
−6
−8

x3

4

4
3
2
1
0
−1
−2
−3
−4

3

x1

2
1
0
−1

6
4

x3

2

−4
−6

4

10

3

8

2

6

1

4

0
−1

2
0

−2

−2

4

4

3

3

2

2

−3

−4

1

x5

x5

0
−2

x4

x4

−2

0

1
0

−1

−1

−3
−3 −2 −1 0

−3
−2 −1

−2

−2
1

x1

2

3

4

−4 −3 −2 −1 0 1 2 3 4

x2

−3 −2 −1

0

x3

1

2

3

−3 −2 −1 0

1

x4

2

3

4

−3 −2 −1 0

1

x5

2

3

4

0

1

x1

2

3

4

−6 −4 −2

0

x2

2

4

(a)

6

−4−3−2−1 0 1 2 3 4 5

x3

−4 −2

0

2

x4

4

6

8

−3 −2 −1 0

1

x5

2

3

4

(a)

FIG. 2. Pair plots for the first 5 dimensions. (a) Data normally distributed as in (37). (b) Data
normally distributed as in (38). We sample 200 points for both cases. We can see that there is a
considerable overlap between the clusters.

improvement of Hartigan’s over Lloyd’s method to be more accentuated in high dimensions.
Thus, we analyze how the algorithms degrade as the number of dimensions increase while
keeping the number of points in each cluster fixed. Consider data from the Gaussian mixture
iid

x ∼ 21 N (µ1 , Σ1 ) + 12 N (µ2 , Σ2 ),
Σ1 = Σ2 = ID ,

µ1 = (0, . . . , 0)> ,
| {z }
×D

µ2 = 0.7(1, . . . , 1, 0, . . . , 0)> .
| {z } | {z }
×10

(37)

×(D−10)

To get some intuition about how separated data points from each class are, we show scatter
plots between the first 5 dimensions in Fig. 2a. Note that the Bayes error is fixed as D
increases, yielding an optimal accuracy of ≈ 0.86. We sample 200 points on each trial. The

results are shown in Fig. 3a. We can see that E H and spectral clustering have practically

the same performance, which is higher than E L -clustering (kernel k-means). Moreover, E H

outperforms k-means and GMM, where the improvement is noticeable specially in high
dimensions. Note that in this setting k-means and GMM are consistent models to the data,
however, energy clustering degrades much less as dimension increases.
Still for a two-class Gaussian mixture, we now allow the diagonal entries of one of the
25

0.95

0.85

0.90

0.75
0.70
0.65
0.60
0.55

0.85

accuracy

accuracy

0.80
EH
EL
spectral
k-means
GMM

0.80
0.75
0.70
0.65
0.60

50

100

150

200

0.55

EH
EL
spectral
k-means
GMM
100

# dimensions

200

300

400

500

600

700

# dimensions

(a)

(b)

FIG. 3. Comparison of E H -clustering, E L -clustering (kernel k-means), spectral clustering, k-means
and GMM in high dimensional Gaussian settings. We plot the mean accuracy versus the number
of dimensions, with error bars indicating standard error from 100 Monte Carlo runs. (a) Data
normally distributed as in (37), with Bayes accuracy ≈ 0.86, over the range D ∈ [10, 200]. (b)
Data normally distributed as in (38), with Bayes accuracy ≈ 0.95, over the range D ∈ [10, 700].

covariances to have different values by choosing
iid

x ∼ 12 N (µ1 , Σ1 ) + 12 N (µ2 , Σ2 ),
µ1 = (0, . . . , 0)> ,
| {z }
×D

µ2 = (1, . . . , 1, 0, . . . , 0)> ,
| {z } | {z }
×10

×(D−10)

Σ1 = ID ,



Σ2 = 

e 10
Σ
0

0
ID−10

e 10 = diag(1.367, 3.175, 3.247, 4.403, 1.249, 1.969, 4.035, 4.237, 2.813, 3.637).
Σ



,

(38)

We simply chose a fixed set of 10 numbers uniformly at random on the interval [1, 5] for the
e 10 , and any other choice would give analogous results. We show pair plots of
diagonal of Σ

this data in Fig. 2b. We sample a total of 200 points from (38) on each trial. The Bayes error

is kept fixed when increasing D yielding an optimal accuracy ≈ 0.95. In Fig. 3b we see that
GMM performs better in low dimensions, but it quickly degenerates as D increases. The
same is true for k-means and E L -clustering. However, E H and spectral clustering remains
much more stable in high dimensions. Notice that a naive implementation of GMM should
not be able to estimate the covariances when D & 100, however, scikit-learn library uses
k-means output as initialization, therefore the output of GMM in this implementation is at
least as good as k-means and the algorithm is more robust in high dimensions.
26

Consider sampling data from the following Gaussian mixture in R20 :
iid

x ∼ 21 N (µ1 , Σ1 ) + 12 N (µ2 , Σ2 ),
µ1 = (0, . . . , 0)> ,
| {z }
×20

µ2 = 12 (1, . . . , 1, 0, . . . , 0)> ,
| {z } | {z }
5

Σ1 = 12 I20 ,

Σ2 = I20 .

(39)

15

The optimal accuracy based on Bayes classification error is ≈ 0.90. We increase the sample
size n ∈ [10, 400] and show the accuracy versus n for the different kernels (34) and (35)

within E H -clustering algorithm, which are compared to k-means and GMM. The results are
in Fig. 4a. Note that for small n all methods are superior than GMM, which slowly catches
up and tend to optimal Bayes, as expected since it is a consistent model to the data. Note
e 1 is as accurate as GMM for large number of points,
also that E H -clustering with kernel K
however, it is superior for small number of points. Still for the same setting, in Fig. 4b we
show the difference in accuracy provided by E H minus E L and E H minus spectral clustering,
e 1 . Note that E H was always superior than kernel k-means and
when using the kernel K
spectral clustering, otherwise there would be points negative values on the y-axis.
Consider the same experiment but now with a lognormal mixture,
iid 1
2

x∼

exp {N (µ1 , Σ1 )} + 12 exp {N (µ2 , Σ2 )} ,

µ1 = (0, . . . , 0)> ,
| {z }
×20

µ2 = 12 (1, . . . , 1, 0, . . . , 0)> ,
| {z } | {z }
5

Σ1 = 12 I20 ,

Σ2 = I20 .

(40)

15

The results are in Fig. 4c. Energy clustering still performs accurately, with any of the
utilized kernels, providing better results than k-means and GMM on this non-normal data.
e 1 still provides the best results for small number of points, but its performance is
The kernel K
eventually achieved by K1/2 , indicating that α ≈ 1/2 in the standard energy distance should

be more appropriate for skewed distributions. In Fig. 4d we show the difference between E H e 1 . Again, the accuracy
clustering to kernel k-means and spectral clustering, with the kernel K
provided by E H is higher than the other methods, although not much higher than spectral
clustering in this example. The two experiments of Fig. 4 illustrate how energy clustering is
more flexible, performing well in different settings with the same kernel, contrary to k-means
and GMM.
In Fig. 5a–c we have complex two dimensional datasets. The two parallel cigars in (a)
have 200 points each. The concentric circles in (b) and (c) have 400 points for each class.
We apply E H -clustering with the kernels (34), (35) and (36). We also consider the best
27

0.90

0.85

0.85

0.80

E H , ρ1

0.75

E H , ρ1/2
E H , ρe1

0.70

50

100

150

200

250

300

350

k-means
GMM

0.70

50

100

150

200

250

# points

# points

(a)

(c)

300

350

400

0.25

difference in accuracy

difference in accuracy

E H , ρe1

0.75

0.55

400

0.25
0.20
0.15
0.10

0.00

E H , ρ1/2

0.80

0.60

0.30

0.05

E H , ρ1

0.65

k-means
GMM

0.65
0.60

accuracy

accuracy

0.90

EH − EL

E H − spectral
50

100

150

200

250

300

350

400

EH − EL

0.20

E H − spectral

0.15
0.10
0.05
0.00

50

100

150

200

250

# points

# points

(b)

(d)

300

350

400

FIG. 4. E H -clustering with kernels (34) and (35) versus k-means and GMM. In both settings Bayes
accuracy is ≈ 0.9. We show average accuracy (error bars are standard error) versus number of
points for 100 Monte Carlo trials. (a,b) Gaussian mixture (39). (c,d) Lognormal mixture (40).
The plots in (c) and (d) consider the difference in accuracy between E H versus E L (kernel k-means)
e1.
and spectral clustering, with the kernel K

(a)

(b)

(c)

(d)

FIG. 5. (a) Parallel cigars. (b) Two concentric circles with noise. (c) Three concentric circles with
noise. (d) MNIST handwritten digits. Clustering results are in Table I and Table II.

28

TABLE I. Clustering data from Fig. 5a–c.
Fig. 5a
ρ1
E H -clustering

spectral-clustering
k-means
GMM

0.705 ± 0.065

ρ1/2 0.952 ± 0.048

ρ1

Fig. 5b

Fig. 5c

0.521 ± 0.005 ρ1

0.393 ± 0.020

ρ1/2 0.522 ± 0.004 ρ1/2 0.486 ± 0.040

ρe2

0.9987 ± 0.0008 ρe1

0.778 ± 0.075 ρe2

0.666 ± 0.007

ρe2

0.557 ± 0.014

ρb1

0.732 ± 0.002 ρb2

0.364 ± 0.004

0.522 ± 0.004 7

0.368 ± 0.005

7

0.595 ± 0.011 7

0.465 ± 0.030

ρb2

0.956 ± 0.020

7

0.550 ± 0.011

7

0.903 ± 0.064

ρb1

7

1.0 ± 0.0

ρb2

0.676 ± 0.002

TABLE II. Clustering MNIST data from Fig. 5d.
{0, 1, . . . , 4}

{0, 1, . . . , 6}

{0, 1, . . . , 8}

{0, 1, . . . , 9}

σ

10.41

10.41

10.37

10.19

ρ1

0.873 ± 0.025

0.731 ± 0.016

0.687 ± 0.016

0.581 ± 0.011

ρ1/2 0.874 ± 0.027

0.722 ± 0.017

0.647 ± 0.017

0.600 ± 0.009

0.847 ± 0.031

0.695 ± 0.023

0.657 ± 0.014

0.584 ± 0.013

Class Subset
parameter

E H -clustering

ρeσ

ρbσ

0.891 ± 0.009 0.759 ± 0.011 0.704 ± 0.011 0.591 ± 0.012
0.769 ± 0.012

0.678 ± 0.014

0.649 ± 0.018

0.565 ± 0.009

k-means

7

0.878 ± 0.010

0.744 ± 0.008

0.695 ± 0.012

0.557 ± 0.012

GMM

7

0.839 ± 0.015

0.694 ± 0.010

0.621 ± 0.009

0.540 ± 0.009

spectral-clustering

ρbσ

kernel choice for each example for spectral clustering. Moreover, we consider k-means and
GMM. We perform 10 Monte Carlo runs for each example. The results are in Table I. For
(a) we initialize all algorithms with k-means++, and for (b) and (c) we initialize at random.
E H has superior performance in every example, and in particular better than the spectral
clustering. In (a) the standard kernel from energy statistics in Euclidean space, K1 and
K1/2 , are able to provide accurate results, however, for the examples in (b) and (c) the
b 1 and K
b 2 provide a significant improvement.
kernel choice is more sensitive, where K

Next, we consider the infamous MNIST handwritten digits as illustrated in Fig. 5d. Each
29

data point is an 8-bit gray scale image forming a 784-dimensional vector corresponding to
the digits {0, 1, . . . , 9}. We compute the parameter

n
1 X
kxi − xj k2 ,
σ = 2
n i,j=1
2

from a separate training set, to be used in the kernels (35) and (36). We consider subsets
of {0, 1, . . . , 9}, sampling 100 points for each class. The results are shown in Table II,

where kernels and parameters are indicated. E H -clustering performs slightly better than
k-means and GMM, however the difference is not considerable. Unsupervised clustering on
MNIST without any feature extraction is not trivial. For instance, the same experiment
was performed in [29] where a low-rank transformation is learned then subsequently used
in subspace clustering, providing very accurate results. It would be interesting to explore
analogous methods for learning a better representation of the data and subsequently apply
E H -clustering.
VIII.

DISCUSSION

We proposed clustering from the perspective of generalized energy statistics, valid for
arbitrary spaces of negative type. Our mathematical formulation of energy clustering reduces to a QCQP in the associated RKHS, as demonstrated in Proposition 2. We showed
that the optimization problem is equivalent to kernel k-means, once the kernel is fixed; see
Proposition 3. Energy statistics, however, fixes a family of standard kernels in Euclidean
space, and more general kernels on spaces of negative type can also be obtained. We also
considered a weighted version of energy statistics, whose clustering formulation establishes
connections with graph partitioning. We proposed the iterative E H -clustering algorithm
based on Hartigan’s method, which was compared to kernel k-means algorithm based on
Lloyd’s heuristic. Both have the same time complexity, however, numerical and theoretical results provide compelling evidence that E H -clustering is more robust with a superior
performance, specially in high dimensions. Furthermore, energy clustering, with standard
kernels from energy statistics, outperformed k-means and GMM on several settings, illustrating the flexibility of the proposed method which is model-free. In many settings, the
iterative E H -clustering also surpassed spectral clustering, which is solves a relaxation of the
original QCQP, and in other settings performed closely but never worse. Note that spec30

accuracy

1.00
0.95
0.90
0.85
0.80
0.75
0.70
0.65
0.60
0.55

EH
EL
spectral
k-means
GMM
0

50

100

150

200

# unbalanced points
FIG. 6. Comparison of energy clustering algorithms to k-means and GMM on unbalanced clusters.
The data is normally distributed as (41), where we vary m ∈ [0, 240], and in each case we do 100
Monte Carlo runs showing the average accuracy with standard error.

tral clustering is more expensive than our iterative method, going up to O(n3 ), and finding
eigenvectors of very large matrices is problematic.
A limitation of the proposed methods for energy clustering is that it cannot handle
accurately highly unbalanced clusters. As an illustration, consider the following Gaussian
mixture:
n1
n1
N (µ1 , Σ1 ) +
N (µ1 , Σ1 ), µ1 = (0, 0, 0, 0)> , µ2 = 1.5 × (1, 1, 0, 0)> ,
2N
2N


1
I2 0
 , n1 = N − m, n2 = N + m, N = 300.
Σ1 = I4 , Σ2 =  2
0 I2
iid

x∼

(41)

We then increase m ∈ [0, 240] making the clusters progressively more unbalanced. We plot
the average accuracy over 100 Monte Carlo runs for each m, with error bars indicating
standard error. The results are shown in Fig. 6. For highly unbalanced clusters we see that
GMM performs better than the other methods, which have basically similar performance.
Based on this experiment, an interesting problem would be to extend E H -clustering algorithm
to account for highly unbalanced clusters.
Moreover, it would be interesting to formally demonstrate cases where energy clustering
is a consistent in the large n limit. A soft version of energy clustering is also an interesting
extension. Finally, kernel methods can benefit from sparsity and fixed-rank approximations
of the Gram matrix, and there is plenty of room to make E H -clustering algorithm more
scalable.
31

ACKNOWLEDGMENTS

We would like to thank Carey Priebe for discussions. We would like to acknowledge
the support of the Transformative Research Award (NIH #R01NS092474) and the Defense
Advanced Research Projects Agencys (DARPA) SIMPLEX program through SPAWAR contract N66001-15-C-4041.

[1] G. J. Székely and M. L. Rizzo. Energy Statistics: A Class of Statistics Based on Distances.
Journal of Statistical Planning and Inference, 143:1249–1272, 2013.
[2] M. L. Rizzo and G. J. Székely. DISCO Analysis: A Nonparametric Extension of Analysis of
Variance. The Annals of Applied Statistics, 4(2):1034–1055, 2010.
[3] G. J. Székely and M. L. Rizzo. Hierarchical Clustering via Joint Between-Within Distances:
Extending Ward’s Minimum Variance Method. Journal of Classification, 22(2):151–183, 2005.
[4] S. Li. k-Groups: A Generalization of k-Means by Energy Distance. PhD Thesis, Bowling
Green State University, 2015.
[5] R. Lyons. Distance Covariance in Metric Spaces. The Annals of Probability, 41(5):3284–3305,
2013.
[6] D. Sejdinovic, B. Sriperumbudur, A. Gretton, and K. Fukumizu. Equivalence of DistanceBased and RKHS-Based Statistic in Hypothesis Testing. The Annals of Statistics, 41(5):2263–
2291, 2013.
[7] S. P. Lloyd. Least Squares Quantization in PCM. IEEE Transactions on Information Theory,
28(2):129–137, 1982.
[8] J. B. MacQueen. Some Methods for Classification and Analysis of Multivariate Observations.
In Proceedings of the 5th Berkeley Symposium on Mathematical Statistics and Probability,
volume 1, pages 281–297. University of California Press, 1967.
[9] E. Forgy. Cluster Analysis of Multivariate Data: Efficiency versus Interpretabiliby of Classification. Biometrics, 21(3):768–769, 1965.
[10] B. Schölkopf, A. J. Smola, and K. R. Müller. Nonlinear Component Analysis as a Kernel
Eigenvalue Problem. Neural Computation, 10:1299–1319, 1998.
[11] M. Girolami. Kernel Based Clustering in Feature Space. Neural Networks, 13(3):780–784,

32

2002.
[12] J. Mercer. Functions of Positive and Negative Type and their Connection with the Theory of
Integral Equations. Proceedings of the Royal Society of London, 209:415–446, 1909.
[13] I. S. Dhillon, Y. Guan, and B. Kulis. Kernel K-means: Spectral Clustering and Normalized
Cuts. In Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge
Discovery and Data Mining, KDD ’04, pages 551–556, New York, NY, USA, 2004. ACM.
[14] I. S. Dhillon, Y. Guan, and B. Kulis. Weighted Graph Cuts without Eigenvectors: A Multilevel
Approach. IEEE Transactions on Pattern Analysis and Machine Intelligence, 29(11):1944–
1957, 2007.
[15] M. Filippone, F. Camastra, F. Masulli, and S. Rovetta. A Survey of Kernel and Spectral
Methods for Clustering. Pattern Recognition, 41:176–190, 2008.
[16] J. A. Hartigan and M. A. Wong. Algorithm AS 136: A k-Means Clustering Algorithm. Journal
of the Royal Statistical Society. Series C (Applied Statistics), 28(1):100–108, 1979.
[17] M. Telgarsky and A. Vattani. Hartigan’s Method: k-Means Clustering without Voronoi.
In Proceedings of the 13th International Conference on Artificial Intelligence and Statistics
(AISTATS), volume 9, pages 313–319. JMLR, 2010.
[18] N. Slonim, E. Aharoni, and K. Crammer. Hartigan’s k-Means versus Lloyd’s k-Means —
Is it Time for a Change? In Proceedings of the 20th International Conference on Artificial
Intelligence, pages 1677–1684. AAI Press, 2013.
[19] N. Aronszajn. Theory of Reproducing Kernels. Transactions of the American Mathematical
Society, 68(3):337–404, 1950.
[20] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola. A Kernel Two-Sample
Test. Journal of Machine Learning Research, 13:723–773, 2012.
[21] C. Berg, J. P. R. Christensen, and P. Ressel. Harmonic Analysis on Semigroups: Theory of
Positive Definite and Related Functions. Graduate Text in Mathematics 100. Springer, New
York, 1984.
[22] J. Shi and J. Malik. Normalized Cust and Image Segmentation. IEEE Transactions on Pattern
Analysis and Machine Intelligence, 22(8):888–905, 2000.
[23] A. Y. Ng, M. I. Jordan, and Y. Weiss. On Spectral Clustering: Analysis and an Algorithm. In
Advances in Neural Information Processing Systems, volume 14, pages 849–856, Cambridge,
MA, 2001. MIT Press.

33

[24] B. Kernighan and S. Lin. An Efficient Heuristic Procedure for Partitioning Graphs. The Bell
System Technical Journal, 49(2):291–307, 1970.
[25] P. Chan, M. Schlag, and J. Zien. Spectral k-Way Ratio Cut Partitioning. IEEE Transactions
on Computer-Aided Design of Integrated Circuits and Systems, 13:1088–1096, 1994.
[26] S. X. Yu and J. Shi. Multiclass Spectral Clustering. In Proceedings Ninth IEEE International
Conference on Computer Vision, volume 1, pages 313–319, 2003.
[27] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel,
P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher,
M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine
Learning Research, 12:2825–2830, 2011.
[28] D. Arthur and S. Vassilvitskii. k-means++: The Advantage of Careful Seeding. In Proceedings
of the Eighteenth annual ACM-SIAM Symposium on Discrete Algorithms, pages 1027–1035,
Philadelphia, PA, USA, 2007. Society for Industrial and Applied Mathematics.
[29] Q. Qui and G. Sapiro. Learning Transformations for Clustering and Classification. Journal
of Machine Learning Research, 16:187–225, 2015.

34

