SoF: Soft-Cluster Matrix Factorization for Probabilistic Clustering
Han Zhao† , Pascal Poupart† , Yongfeng Zhang§ and Martin Lysy‡†

David R. Cheriton School of Computer Science, University of Waterloo, Canada
Department of Computer Science and Technology, Tsinghua University, China‡
Department of Statistics and Actuarial Science, University of Waterloo, Canada†‡{han.zhao, ppoupart, mlysy}@uwaterloo.ca, § zhangyf07@gmail.com§

Abstract
We propose SoF (Soft-cluster matrix Factorization), a probabilistic clustering algorithm which softly assigns each data
point into clusters. Unlike model-based clustering algorithms,
SoF does not make assumptions about the data density distribution. Instead, we take an axiomatic approach to define 4
properties that the probability of co-clustered pairs of points
should satisfy. Based on the properties, SoF utilizes a distance measure between pairs of points to induce the conditional co-cluster probabilities. The objective function in
our framework establishes an important connection between
probabilistic clustering and constrained symmetric Nonnegative Matrix Factorization (NMF), hence providing a theoretical interpretation for NMF-based clustering algorithms.
To optimize the objective, we derive a sequential minimization algorithm using a penalty method. Experimental resultson both synthetic and real-world datasets show that SoF significantly outperforms previous NMF-based algorithms andthat it is able to detect non-convex patterns as well as cluster
boundaries.

Introduction
Clustering algorithms have played an important role in various research tasks, including but not limited to document
clustering (Steinbach et al. 2000; Xu, Liu, and Gong 2003),
image segmentation (Shi and Malik 2000) and gene clustering (Thalamuthu et al. 2006). Recently, probabilistic clustering has gained in popularity and application (He et al. 2011;
Jain 2010; Aggarwal and Reddy 2013) because the probabilistic nature of the clusters facilitates their integration into
several tasks including latent semantic indexing (Ding, Li,
and Peng 2008), probabilistic pre-training in deep learningtechniques (Erhan et al. 2010), etc. Most of the existing
probabilistic clustering algorithms attempt to infer the posterior cluster distribution by assuming a family of density
distributions, typical examples include the Gaussian mixture
model, the Dirichlet process (Biernacki, Celeux, and Govaert 2000), the Wishart-Dirichlet process (Vogt et al. 2010),
etc. However, the assumed family may not include the distribution of the data, hence a flexible and data-driven approachis desirable.
Copyright c 2015, Association for the Advancement of Artificial
Intelligence (www.aaai.org). All rights reserved.

In this paper, we propose a novel probabilistic clustering
algorithm based on Soft-cluster Matrix Factorization (SoF).
Different from previous approaches, SoF does not make assumptions about the data density distribution. Instead, SoF
employs a distance function between pairs of input points
to induce their co-cluster probabilities. By directly considering the co-cluster probabilities, SoF can model clusters
with various patterns, including both convex and non-convex
clusters. The input to SoF consists of a set of data points, a
distance function L over these data points, and the numberof clusters, K. As we will see later, the choice of distancefunction L reflects our prior knowledge about the clusteringtask.
SoF has a close connection with constrained symmetric NMF. In the past few years, algorithms for solving
NMF have been studied extensively (Seung and Lee 2001;
Kim and Park 2008a; 2008b; 2008c; Kuang, Park, and Ding
2012) and applied to clustering tasks. However, to the bestof our knowledge, the connection between NMF and probabilistic clusterings is not well understood. The objectivefunction in our framework reveals an intrinsic connection
between probabilistic clustering and constrained symmetric
NMF, hence giving a theoretical justification for the application of symmetric NMF techniques to clustering tasks. Furthermore, we will show that the true pairwise co-cluster matrix for probabilistic clustering is actually a completely positive matrix, which is NP-hard to determine (Berman andShaked-Monderer 2003). Our clustering algorithm can thenbe viewed as approximating the completely positive matrix
in a relaxed cone space, which is a better approximation thanother existing NMF based algorithms. To solve SoF efficiently, we derive a sequential minimization framework to
optimize the objective. Experiments on both synthetic andreal-world datasets show that SoF is able to find non-convex
clusters due to the absence of any density assumption andits probabilistic nature can be used to detect boundaries between multiple clusters.

Preliminaries
Given a set of input points I = {vi ∈ Rd | i = 1, . . . , N }and the number of clusters K, the task of probabilistic clustering is to associate a probability vector p ∈ RK to eachpoint v ∈ I, where the jth component of p represents theprobability that v is assigned to the jth cluster. The proba-bility vector p should satisfy p  0 and 1T p = 1.1
Let V be a continuous random vector in Rd for input
points and C be a categorical random variable which takes
a cluster label from {1, . . . , K}. We use p(V, C) to denote the joint probability distribution over domain(V) ×domain(C). Each sample from the joint probability distribution is a pair (v, c), where v is the point we observe and cdenotes the cluster label associated with v. We assume that
all the pairs (point and its label) are sampled from a mixture model independently. More precisely, we assume there
are K components where each component corresponds to acluster. Each component can generate observed data basedon its conditional distribution p(V|C = k), k = 1, . . . , K
independently. Unlike traditional mixture models, we do
not assume any fixed form for the conditional distributionsp(V|C = k). Moreover, for different k’s, the conditional
distribution can take different forms.

Soft-cluster Matrix Factorization
We formally define a probabilistic clustering function as follows:
Definition 1. Given an input set I = {vi ∈ Rd | i =1, . . . , N } with some fixed dimension d, a distance function L : Rd × Rd → R+ that is reflexive, symmetric andpositive definite2 , and a fixed number of clusters K, a probabilistic clustering algorithm is a function that takes I, L, K
as inputs and for each v ∈ I, outputs a vector pv ∈ SK−1 ,where SK−1 represents the K − 1 dimensional probability
simplex.
Given I, L, K, a probabilistic clustering algorithm
should, for each v ∈ I, output a probability vector which
is close to the true posterior probability vector, i.e., pv ≈(p(C = 1|V = v), . . . , p(C = K|V = v))T . The choice ofdistance function L(·, ·) reflects our prior knowledge about
the clustering task, where the basic assumption is that points
which are close in terms of L tend to reside in the same cluster and points which are far away from each other in termsof L tend to reside in different clusters.

1. Boundary property: ∀v1 , v2 ∈ Rd , if L(v1 , v2 ) = 0, thenpC (v1 , v2 ) = 1; if L(v1 , v2 ) → ∞, then pC (v1 , v2 ) →0.
2. Symmetry property: ∀v1 , v2 ∈ Rd , pC (v1 , v2 ) =pC (v2 , v1 ).
3. Monotonicity property: ∀v1 , v2 , v3 ∈ Rd , pC (v1 , v2 ) ≤pC (v1 , v3 ) if and only if L(v1 , v2 ) ≥ L(v1 , v3 ).
4. Tail property: Given ∀v1 , v2 , v3 ∈ Rd ,∂pC (v1 ,v3 )∂L(v1 ,v3 )∂pC (v1 ,v2 )∂L(v1 ,v2 )≤if and only if L(v1 , v2 ) ≥ L(v1 , v3 )The boundary property forces pC (·, ·) to be a valid probability metric and the symmetry property allows the argumentsof pC (·, ·) to be interchangeable. The monotonicity property
formalizes our intuition that closer points tend to be more
likely to reside in the same cluster. This property also suggests that L(·, ·) should be well-chosen to reflect our priorknowledge about the data. The tail property can be interpreted as follows: changes in the distance of two points will
have an impact on the co-cluster probability that diminisheswith the increase of the distance between them.
Proposition 1. Given a distance function L : Rd × Rd →R+ , there exists a family of co-cluster probability functionspC (v1 , v2 ) = e−cL(v1 ,v2 ) which satisfy the boundary, symmetry, monotonicity and tail properties simultaneously for
any constant c > 0.
Proof. The boundary property is satisfied since e−c0 = 1
and e−c∞ = 0. The symmetry property follows from the factthat L is symmetric. The monotonicity property holds sincee−cL(v1 ,v2 ) ≤ e−cL(v1 ,v3 ) ⇐⇒ L(v1 , v2 ) ≥ L(v1 , v3 ).
The tail property can be derived as follows:∂pC (v1 , v2 )∂pC (v1 , v3 )≤∂L(v1 , v2 )∂L(v1 , v3 )⇐⇒ −ce−cL(v1 ,v2 ) ≤ −ce−cL(v1 ,v3 )⇐⇒ e−cL(v1 ,v2 ) ≤ e−cL(v1 ,v3 )⇐⇒ L(v1 , v2 ) ≥ L(v1 , v3 )Co-cluster Probability
The distance function L : Rd × Rd → R+ measures how
close two points are. One of the key questions in clustering is
how can we obtain a similarity measure between input points
in terms of the distance function L? More specifically, whatis the probability that two points belong to the same cluster?
We propose the following four basic properties to find appropriate expressions for the probability that a pair of points
is co-clustered based on the distance function L. Hence,∀v1 , v2 ∈ Rd , we will use pC (v1 , v2 ) to denote the empirical probability that v1 and v2 reside in the same cluster.
A co-cluster probability pC (·, ·) should satisfy the following
four properties:
1
We use  to denote the element-wise greater than between
matrices of the same dimension.
2
The distance function L(·, ·) is not required to be a metric.

For the specific choice where L is the squared Euclidean
distance, our construction of the co-cluster probability function coincides with the Gaussian kernel. Similarly, when L
is the Euclidean distance, the co-cluster probability function
is equivalent to the Laplacian kernel. Note however that ourframework is more general in the sense that one may use
other kinds of distance functions L that may not lead to wellknown kernels. In e−cL(v1 ,v2 ) we have a scaling parameterc. When c → 0, the probability that any two points belong to
the same cluster tends to 1 and a single global cluster holds
all the data points; when c → ∞, the co-cluster probabilityof any two points tends to 0 and each point forms its own
cluster. Hence, the clustering quality is very sensitive to the
choice of scaling parameter c. Later we will discuss a strategy to avoid manually tunning the scaling parameter c whileat the same time achieving the scaling invariant property ofclustering (Kleinberg 2003). Henceforth, for simplicity of
the notation, we will implicitly assume c = 1.

Probabilistic Clustering
We use P ∈ RN ×N to denote the soft-cluster matrix, where
Pij , pC (vi , vj ) = e−L(vi ,vj ) is the empirical probability
that two points are clustered together. In our setting, P acts
as the refined prior knowledge extracted from the raw input
to guide the probabilistic clustering algorithm in the learning
process.
Based on our probabilistic model, given vi , vj , i 6= j,
the true probability that vi , vj belong to the same cluster,Pr(vi ∼ vj ), is given by
Pr(vi ∼ vj ) =KXp(ci = cj = k|vi , vj )k=1
K
ζ1 X p(vi , ci = k)p(vj , cj = k)=p(vi , vj )

k=1
K
ζ2 X

=

p(ci = k|vi )p(cj = k|vj ) = pTvi pvj

k=1

where ζ1 and ζ2 follow from the i.i.d. assumption of the data
generation process.
Given the joint probability distribution p(V, C), we cancompute the conditional probability distribution of cluster label C given V by p(C|V) = p(V,C)
p(V) . Let W =×K[Wij ]N ×K ∈ RN
, where Wij = p(ci = j|vi ). We+
rewrite Pr(vi ∼ vj ) as:Pr(vi ∼ vj ) =K
X

Wik Wjk = Wi: W:jT

(1)

k=1

where Wi: represents the ith row vector of matrix W . Let’s
define Pr(V ∼ V0 ) , [Pr(vi ∼ vj )]N ×N and express it inmatrix form:
Pr(V ∼ V0 ) = W W T(2)The ith row of W corresponds to the conditional probability vector that vi belongs to different clusters, i.e.,
Wi: = pTvi . We see that the true co-cluster probability matrix Pr(V ∼ V0 ) is both nonnegative and symmetric positive semidefinite. Furthermore, it can be factorized as the
product of a nonnegative matrix and its transpose. Such amatrix is known as the Completely Positive matrix (C.P. matrix) in the literature (Berman and Shaked-Monderer 2003;
Dickinson and Gijben 2014). C.P. matrices are containedin the intersection of the positive semidefinite cone andthe nonnegative cone, but they are a strict subset of thisclass (Dickinson and Gijben 2014). We will use the C.P.
cone to denote the set of C.P. matrices. It has recently been
proven that determining whether a given matrix is in the C.P.
cone is NP-hard (Dickinson and Gijben 2014). From thisperspective, our construction of the soft-cluster matrix is actually a relaxation that allows us to search in the intersectionof the nonnegative cone and the positive semidefinite cone,
which is a strict superset of the C.P. cone.

Our goal is to learn the matrix W from the inputs I, L, K.
Given I, L, K, the soft-cluster matrix P , which satisfies thefour properties and admits the form e−L(·,·) , is a relaxation
and approximation of the true co-cluster matrix Pr(V ∼V0 ). Therefore we obtain a natural learning paradigm to estimate the optimal W , denoted by W ∗ :minimizeW
subject to||P − W W T ||2F
N ×K
W ∈ R+, W 1K = 1N

(3)

where 1K and 1N are K and N dimensional column vectors
whose components are all 1; || · ||F represents
the FrobeniusPnorm of a matrix where ||A||F = ( i,j A2ij )1/2 . Before
tackling the optimization problem above, we discuss someof its properties and its relations to other clustering algorithms.
1) SoF is closely related to other clustering algorithms, including Kernel Kmeans (Dhillon, Guan, and Kulis
2004), Spectral Clustering (Ding, He, and Simon 2005)
and Symmetric Nonnegative Matrix Factorization (SymNMF) (Kuang, Park, and Ding 2012), but it also exhibits significant differences. We list the connections and differencesamong those methods in Table 1. In summary, the matrix
P used in SoF is nonnegative as well as symmetric positive
semidefinite (S.P.D.), which distinguishes it from all otherthree methods. Meanwhile, the linear constraints W 1 = 1
together with the nonnegativity constraint in SoF make thesolution clearly interpretable as a probabilistic assignmentof points to clusters, while all other three methods seek tofind data representations in either kernel spaces or reduced
eigenspaces.
2) The decision version of the optimization problem in (3)is closely connected to the strong membership problem for
the C.P. cone, which asks given an S.P.D matrix P which isnonnegative, whether there exists a nonnegative matrix W
such that P = W W T . It’s then easily seen that the decision version of (3) is essentially a restricted version of thestrong membership problem for the C.P. cone, which is NPhard. We conjecture that the decision version of (3) is also
NP-hard because of its close connection to the strong membership problem for the C.P. cone. If this is the case, the optimization problem formulated in (3) is also NP-hard. However, the validation of this conjecture is out of the scope of
this paper and we leave it for future work.
3) The optimization problem (3) is non-convex in generalbecause the objective function contains a fourth order polynomial in W whose Hessian matrix cannot be guaranteed
to be positive semi-definite. The nonnegative constraint together with the N linear constraints make it even harder tosolve exactly because it is known that solving NMF exactly
is NP-hard (Vavasis 2009). Therefore, instead of looking fora global optimal solution, we seek to find a local optimal solution that is tractable to compute and at the same time close
to a global optimum.
4) The optimal solution W ∗ is not unique. Given an optimal solution W ∗ , we can construct another optimal solutionŴ ∗ = W ∗ Π, where Π is a K × K permutation matrix.3 The
3A permutation matrix is a square binary matrix that has exactly
one entry 1 in each row and each column and 0s elsewhere.

effect of right multiplication of Π is to permute the labels of
the clusters. Considering all possible permutations of cluster
labels, there are at least K! optimal solutions. However, any
one of the optimal solutions suffices since we do not careabout the exact labeling of each cluster.
5) In clustering, we often prefer a sparse solution matrix
W ∗ since only the data points near the boundary of multiple
clusters may have a dense probability vector, which meansthat the majority of the row vectors in W ∗ should be sparse.
This observation will help us to construct an equivalent characterization of (3) in the next section. Furthermore, we canchoose a sparse initial matrix W (0) as a starting point for ouriterative algorithm introduced in the next section to optimize
the objective function based on the sparsity nature of W .

To apply the sequential minimization framework, we begin by deriving the gradient for the unconstrained objectivefunction. Throughout this paper we will use the numerator
layout (Minka 2000) to compute the differentials, then thegradient is simply the transpose of the corresponding differential. First consider the differential of ||P − W W T ||2F withrespect to W :d||P − W W T ||2F = dtr{(P T − W W T )(P − W W T )}= 2tr{(W W T − P T )dW W T }= 4tr{W T (W W T − P T )dW }Note that we are using the fact tr(AB) = tr(BA) together
with tr(AT ) = tr(A). Hence the gradient of the first term in(4) is:Optimization Algorithm
In this section, we focus on solving the problem in (3) formulated as a constrained norm-minimization problem. As
discussed in the previous section, it is difficult to solve
this problem exactly. Hence, we seek an approximation
procedure by converting the original constrained optimization problem into a sequence of unconstrained optimization
problems by using a penalty method. More concretely, wetransform the two constraints in (3) into penalty functions
that are added to the objective to penalize solutions that are
not feasible:P
minimizeW ||P − W W T ||2F − λ1 ij min{0, Wij }+λ2 ||W 1K − 1N ||22(4)
The second term in (4) is a hinge-loss that penalizes elements in W which do not lie in the nonnegative
cone. OneP
can also use the log-barrier function −λ1 ij log(Wij ) as
the penalty term. Here we choose the hinge-loss rather thanthe log-barrier because we seek a sparse solution (the logbarrier function often pushes the solution away from theboundary of the feasible region until the coefficient λ1 → 0).
The third term in (4) is an `2 regularization term to penalize
rows of W that violate the linear constraints.
In (4), λ1 , λ2 ∈ R+ determine the tradeoff between theaccuracy of the approximation procedure and the feasibilityof the solution. With λ1 and λ2 → ∞, any violations of theconstraints will be greatly penalized and the solution is guaranteed to be in the feasible region. On the other hand, theobjective term ||P − W W T ||2F may be dominated by thesetwo penalty terms, causing the solution to be inaccurate.
To combine the best of two worlds, we propose a sequential minimization framework: for each fixed pair (λ1 , λ2 ),we optimize (4) to find a local minimum solution W (λ1 , λ2 )which depends on the current (λ1 , λ2 ), then we increase theparameters (λ1 , λ2 ) by a factor µ > 1 and restart the optimization with W (λ1 , λ2 ) as the initial point. To increase
the approximation accuracy when λ1 and λ2 get larger, weincrease the number of steps per iteration, ensuring that thepenalty terms will vanish as the solution becomes feasible,
and the objective term will get fully optimized. The wholeprocess is repeated until the current solution is a stationary
point.

∇W {||P − W W T ||2F } = 4(W W T − P )W

(5)

Using the same trick we can compute the gradients of thesecond term and third term in (4) as:X∇W {−λ1
min{0, Wij }} = −λ1 (IW ≺0 1N 1TK ) (6)ijand∇W {λ2 ||W 1K − 1N ||22 } = 2λ2 (W 1K 1TK − 1N 1TK ) (7)where is the Hadamard product of two matrices, i.e., theelement-wise product, and IW ≺0 is the indicator function
Ix<0 applies to every element of W .
The gradient of the unconstrained optimization objectiveis the combination of the three gradients above. Once we obtain a closed form solution to compute the gradient matrix,
we can apply various kinds of unconstrained minimization
algorithms, e.g., gradient descent, conjugate gradient or LBFGS (Andrew and Gao 2007) to obtain a local minimumof (4). For illustrative purposes, we provide an algorithm using gradient descent in Alg. 1. Assuming the algorithm will
perform τ iterations in each inner loop, the time complexityof Alg. 1 is O(τ N 2 K logµ (λ/)).

Experiments
We first review how to construct a reasonable distance function L required in Alg. 1. To demonstrate our algorithm, wefirst apply SoF on synthetic data, which consists of both convex and non-convex patterns, and then we compare SoF withmultiple clustering algorithms on real-world datasets.

Distance Function
Our probabilistic clustering algorithm requires a distancefunction L(·, ·) over all pairs of points as an input. As mentioned before, the choice of scale parameter c will have an
impact on the final clustering and furthermore is problem
dependent. To avoid an annoying tuning process, we propose to use a relative distance L0 (·, ·) induced from the absolute distance L(·, ·). More specifically, for any two points
L(v ,v )vi , vj ∈ I, we define L0 (vi , vj ) as L0 (vi , vj ) = √σini√σjjn ,where σin is the distance between vi and the nth closest
point of vi measured by L. Intuitively, σin defines the scopeTable 1: Connections and differences among clustering algorithms

Objective
Property
Constraint

Kernel Kmeansmin ||K − W W T ||2F
K is S.P.D.
W T W = I, W  0Spectral Clusteringmin ||L − W W T ||2F
L is graph Laplacian, S.P.D.
WTW = I

Algorithm 1 Sequential Minimization to Unconstrained
Problem
Input: I = {v1 , . . . , vN }, distance function L, # of clusters(0)K, scale factor c > 0, initial penalty parameters λ1 ,(0)
λ2 > 0, step-factor µ > 1, initial learning rate 0 <(0)
γ < 1, threshold 0 < 1 , 2  1.
Output: K dimensional probability vector pv ∀v ∈ I.
1: Build distance matrix D such that Dij = L(vi , vj ).
2: Build soft-cluster matrix P such that Pij = e−cDij .
×K
3: Initialize matrix W (0) ∈ RN.
+
4: repeat
(t)
(t−1)
5:
W0 ← W∗6:
repeat(t)
7:Compute the gradient ∇W f (Wl ) according to(5), (6) and (7).
(t)
(t)
(t)
8:
Wl+1 ← Wl − γ (t) ∇W f (Wl ).
9:(t)

SymNMFmin ||A − W W T ||2F
A is similarity matrix
W 0

SoF
min ||P − W W T ||2F
P is nonnegative, S.P.D.
W  0, W 1 = 1

Perona 2004) and generated by (Pei and Zaı̈ane 2006)4 .
Due to the probabilistic nature of SoF, for each point v,
we pick its label by maximizing the conditional probabilityp(c = k|v), c = 1, . . . , K. We break ties randomly.
Clustering Quality The 6 patterns shown in Fig. 1 havedifferent cluster structures, and some of them are not convex, e.g., 1, 2, 3 and 6. The clustering results are shown in
Fig. 1, where different colors correspond to different clusters
found by SoF. Since SoF does not make any density assumption, it is able to detect various cluster structures, including
convex and non-convex patterns. The fifth pattern in Fig. 1
is generated from 2 Gaussians. However, with the numberof clusters K set to be 3, our algorithm is still able to find areasonable clustering and to quantify the uncertainties along
the cluster gap (as shown in Fig. 2).

(t)until ||Wl+1 − Wl || < 1 .
(t)

(t)

W∗ ← W L
(t+1)
(t)
11:
λ1
← λ1 × µ.
(t+1)
(t)
12:
λ2
← λ2 × µ.
(t+1)
(t)
13:
γ
← γ /µ.
(t)(t)
14: until λ1 > 1/2 and λ2 > 1/2
15: return W ∗10:of vi , which, to some extent, converts the absolute distance
between two points into a relative distance. This transformation will also give us the scaling invariant property in thesense that enlarging or shrinking distances between pairs of
points by a constant factor will not affect the clustering result. More specifically, consider the distance L(·, ·) between
a pair of points that have been scaled by a constant factor
α > 0, the relative distance will still remain the same:

Figure 1: 6 synthetic patterns, where the # of clusters K fromleft to right and top to bottom is (1) K = 3, (2) K = 5, (3)K = 3, (4) K = 3, (5) K = 3, (6) K = 2.

15

15

0.03

5

αL(vi , vj )L(vi , vj )L (vi , vj ) = √=√ √= L0 (vi , vj )√ασin ασjn
σin σjn
0

15

0.035
10

0.025
0.02

0

0.8
10

0.7
0.6

5
0

0.4
0.3

−5
0.01

−10

0.005

−15
−15

−10

−5

0

5

10

15

0

−10
−15
−15

0.1
−10

−5

0

5

10

15

1.4

In our experiments, we choose L(·, ·) to be the Euclidean
distance between each pair of points and set n = 10. Different values of n will still affect the clustering, but it will notaffect it too much in practice.

10

1.2
1

5

0.8

0

0.6
0.5
0.4
−5

0.3
0.2

−10

0.1

0

−15
−15

1.5

15

−10

−5

0

5

10

15

10

10

5

1

0.8
0.7

5

0.6

0

0

−5

−5

0.5
0.4

0.5

0.3

0.4
−10

0.2
−10

Synthetic Dataset

−5

0

5

10

15

0

−10

0.2

−10

0.1
−15
−15

−10

−5

0

5

10

15

0

−15
−15

−10

−5

0

5

10

Figure 2: Entropy graphs for 6 patterns in Fig. 1.
4

0

0.9

0.6
−5

−15
−15

To visualize the performance of SoF, we perform experiments on the synthetic data set used in (Zelnik-Manor and

0.7

0

0.2

15

15

0.8

5

0.5

0.015
−5

0.9
10

http://webdocs.cs.ualberta.ca/˜yaling/
Cluster/Php/data_gen.php

150Clustering Boundary by Quantifying Uncertainties
Given a probability vector p ∈ SK−1 , we use entropy tomeasure the uncertainty of cluster assignments of the corPK
responding data points: entropy(p) = − i=1 pi log pi .
Hence, points on the cluster boundaries will have a higher
entropy while those that are inside a cluster should havelower entropy. Fig. 2 shows the entropy graphs wherebrighter colors, including yellow and red, indicate a higher
entropy value while dark colors such as blue and dark green
indicate a lower entropy value. We can visualize the cluster
boundaries by checking points associated with brighter colors in convex patterns, e.g. 4 and 5. Furthermore, it can beseen that most of the points (with dark blue colors) belong
to a single cluster, indicating that the solutions obtained by
SoF are sparse, which is a desirable property in clustering.

Real-world Dataset
Table 2: Statistics of Datasets
Data setsblood (BL)breast t (BT)glass (GL)iris (IRIS)
ecoli (ECO)
satimage (IMG)
pendigits (DIG)# Objects
748
106
214
150
336
4435
10992

# Attributes
4
9
9
4
7
36
16

# Classes
2
6
6
3
8
6
10

Experimental Setup We use 7 real-world data sets fromthe UCI Machine Learning Repository5 . The statistics of thedatasets are summarized in Table 2. More detailed information about these data sets can be found at the UCI Machine
Learning Repository. The points in each dataset are associated with true labels, which will be used as ground truth forthe evaluation of different clustering algorithms.
To evaluate the performance of SoF and compare it withother hard clustering algorithms, we assign each point to thecluster with the highest probability: lv = arg maxc p(c|v).
It is commonly admitted that the validation of clustering
structures is the most difficult task of clustering analysis(Jain and Dubes 1988). In our experiments we use three robust measures to evaluate clustering performance, namely
purity, rand index and accuracy (Jain and Dubes 1988;
Cai, He, and Han 2011; Murphy 2012).
1. Purity. Let Nij be the number of objects in cluster i thatPC
belong to class j, and let Ni =j=1 Nij be the total
number of objects in cluster i. Define pij = Nij /Ni ;this is the empirical probability that a point in cluster
i is from class j. The purity of a cluster is defined aspi , maxj pij , and the overall purity of a clustering isdefined as
X Nipipurity ,Ni
The purity ranges from 0 to 1, and the higher the purity,
the better the clustering result.
5archive.ics.uci.edu/ml/datasets.html

2. Rand Index. Let U = {u1 , . . . , uR } and V ={v1 , . . . , vC } be two different partitions of the N data
points, i.e., two different flat clusterings. We assume V tobe the ground truth clustering, which is induced by classlabels. Now define a 2 × 2 contingency table, containing the following numbers: TP, TN, FP and FN. TP is thenumber of pairs that are in the same cluster in both U andV (true positives); TN is the number of pairs that are in
the different clusters in both U and V (true negatives); FN
is the number of pairs that are in the different clusters in
U but the same cluster in V (false negative); and FP is
the number of pairs that are in the same cluster in U but
different clusters in V (false positive). The Rand Index isdefined as:
TP+TN
R,
TP+FP+FN+TN
Rand index can be interpreted as the fraction of correctly
clustered pairs. Clearly 0 ≤ R ≤ 1, and the larger therand index, the better the clustering.
3. Accuracy. Let N be the total number of objects. For the
ith object, we use ui to denote the label generated by
the clustering algorithm and vi to denote the ground truthclass label associated with the ith object. The accuracy isdefined asδ(vi , map(ui ))accuracy =N
where δ(·, ·) is the Kronecker delta function, map(·) isthe permutation that maps each cluster label generated by
the clustering algorithm to a class label of the data corpus. The best mapping can be obtained by Kuhn-Munres’algorithm (Lovász and Plummer 1986). The higher the accuracy, the better the clustering result.
We compare the clustering qualities based on purity, rand
index and accuracy of the following clustering algorithms:
1. Standard K-means (K-means) and Spherical K-means (Smeans). For two data points vi and vj , standard K-means
uses the Euclidean distance ||vi − vj ||2 as its distancefunction while spherical K-means uses 1 − cos(vi , vj ) as
its distance function.
2. Gaussian Mixture Model (GMM). EM algorithm is used
to fit data density distributions and the number of components is given by the true number of clusters in the dataset.
3. Spectral clustering-Normalized Cut (SP-NC) and Spectral clustering-Ng (SP-NJW). There are two streams ofspectral clustering algorithms, namely normalized cut,proposed by Shi and Malik, and Ng-Jordan algorithm,
proposed by Ng et al.. In both algorithms, K-means isused as the final procedure to obtain the low-dimensional
clustering.
4. Standard NMF and Regularized NMF (RNMF). Nonnegative matrix factorization based methods, including
NMF and its regularized variant, work on a sample matrix directly to obtain a low dimensional representationof the original high dimensional input data. These algorithms interpret each data point by its low dimensional approximation and assign each data point its cluster label bymaximal approximation coefficient in the lower dimensional space. In our experiments, we use the alternating
non-negative least square algorithm (Kim and Park 2008a;
2008c). Based on standard NMF, regularized NMF adds
two Frobenius regularized terms to penalize non-negative
matrix factorization with large elements in the two resulting matrices.
5. GNMF. Cai et al. proposed GNMF to take neighboring
relationships into account by adding a graph-theoretic
penalty term to standard NMF.
6. SymNMF. Kuang, Park, and Ding proposed SymNMF asa Newton-like algorithm to solve symmetric non-negative
matrix factorization problems, including graph clustering
problems. Different from our probabilistic clustering algorithm, SymNMF factorizes the Graph Laplacian Matrix (Shi and Malik 2000; Ng et al. 2002) and interprets
the result as NMF does.
7. SoF. Our proposed probabilistic clustering algorithm using sequential gradient descent minimization.
Table 3: Average purity over 20 runs on 7 UCI data sets.
K-means
S-means
GMM
SP-NC
SP-NJW
NMF
RNMF
GNMF
SymNMF
SoF

BL
0.76
0.76
0.77
0.78†
0.76
0.76
0.76
0.76
0.76
0.76

BT
0.43
0.43
0.46
0.42
0.40
0.43
0.28
0.44
0.47
0.51†

GL
0.56
0.57
0.52
0.38
0.38
0.58
0.52
0.55
0.60
0.64†

IRIS
0.87
0.92
0.89
0.86
0.59
0.79
0.84
0.85
0.88
0.95

ECO
0.80
0.76
0.81
0.81
0.80
0.73
0.70
0.77
0.81
0.85†

IMG
0.70
0.63
0.74
0.30
0.48
0.57
0.56
0.71
0.77†
0.75

DIG
0.71
0.71
0.70
0.11
0.13
0.47
0.45
0.71
0.76
0.82†

Table 5: Average accuracy over 20 runs on 7 UCI data sets.
K-means
S-means
GMM
SP-NC
SP-NJW
NMF
RNMF
GNMF
SymNMF
SoF

BL
0.73
0.61
0.66
0.73
0.56
0.68
0.71
0.72
0.76
0.76

BT
0.34
0.39
0.43
0.41
0.42
0.39
0.27
0.36
0.42
0.48†

GL
0.51
0.52
0.46
0.37
0.33
0.50
0.51
0.44
0.46
0.47

IRIS
0.86
0.90
0.86
0.85
0.56
0.79
0.84
0.81
0.87
0.94

ECO
0.60
0.56
0.68
0.61
0.58
0.66
0.59
0.56
0.66
0.74†

IMG
0.63
0.61
0.69
0.28
0.37
0.53
0.52
0.54
0.68
0.71†

DIG
0.68
0.67
0.67
0.11
0.12
0.44
0.41
0.67
0.67
0.82†

ous NMF based clustering methods, including NMF, RNMF,
GNMF and SymNMF because of SoF’s sequential minimization framework. One possible direction to make SoF
more scalable is to utilize the potential sparse structure of
the empirical pairwise co-cluster matrix.
Since the datasets from different domains present different patterns, our probabilistic clustering algorithm is not biased to a specific domain and it achieves superior resultson all three measures in nearly all of the clustering tasks.
For some tasks, the comparative algorithms may happen to
make the right assumption about the underlying data density
and thus yield better results (e.g., spectral clustering is thebest for blood), but these assumptions may hurt the performance when they do not hold in other datasets (e.g., spectralclustering performs badly on satimage and pendigits). Furthermore, SoF consistently outperforms other NMF basedalgorithms, including NMF, RNMF and GNMF, on almost
all the datasets.

Conclusions
Table 4: Average rand index over 20 runs on 7 UCI data sets.
K-means
S-means
GMM
SP-NC
SP-NJW
NMF
RNMF
GNMF
SymNMF
SoF

BL
0.60
0.52
0.57
0.64
0.51
0.56
0.59
0.60
0.62
0.64

BT
0.71
0.76
0.76
0.74
0.71
0.65
0.38
0.68
0.77
0.79†

GL
0.69
0.69
0.64
0.43
0.53
0.71
0.65
0.70
0.70
0.73†

IRIS
0.86
0.92
0.88
0.85
0.60
0.81
0.85
0.85
0.88
0.93

ECO
0.81
0.78
0.85
0.81
0.80
0.81
0.78
0.80
0.83
0.85

IMG
0.81
0.81
0.83
0.32
0.74
0.80
0.79
0.77
0.83
0.86†

DIG
0.91
0.91
0.91
0.11
0.41
0.86
0.86
0.91
0.91
0.94†

Results and Analysis Experimental results on 7 UCI data
sets are reported in Tables 3, 4 and 5. We run each clustering algorithm 20 times on each data set and report the average performance. The algorithms which achieve the highest mean over 20 runs are highlighted with bold font andthose that are statistically better according to the Wilcoxon
signed-rank test (Wilcoxon and Wilcox 1964) with p-value≤ 0.05 are marked with a dagger, †. Among all the algorithms, K-means and S-means are the fastest. As the algorithm is currently implemented, SoF is slower than previ-

In this paper we propose a novel probabilistic clustering
algorithm whose objective is axiomatically derived from a
set of properties characterizing the co-cluster probability interms of pairwise distances. Our approach can be viewed
as a relaxation of C.P. programming, which intrinsically reveals a close relationship between probabilistic clusteringand symmetric NMF-based algorithms. We further design asequential minimization framework to find a local minimum
solution. Experiments on synthetic and real-world datasetsshow encouraging results. Future work includes speeding up
our sequential minimization framework by utilizing the potential sparsity structure of the empirical pairwise co-cluster
matrix so that it scales linearly with the non-zero elementsin the matrix.

Acknowledgments
The authors are grateful to Da Kuang for helpful discussions
as well as Hongfu Liu for their datasets. The work done
by Han Zhao and Pascal Poupart was supported by NSERC
Discovery and Engage grants. Han Zhao was also supported
by David R. Cheriton Scholarship. Yongfeng Zhang was
supported by Baidu and IBM PhD Fellowship program.

References
Aggarwal, C. C., and Reddy, C. K. 2013. Data Clustering:
Algorithms and Applications. CRC Press.
Andrew, G., and Gao, J. 2007. Scalable training of l 1regularized log-linear models. In Proceedings of the 24th international conference on Machine learning, 33–40. ACM.
Berman, A., and Shaked-Monderer, N. 2003. Completely
positive matrices. World Scientific.
Biernacki, C.; Celeux, G.; and Govaert, G. 2000. Assessing
a mixture model for clustering with the integrated completed
likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence 22(7):719–725.
Cai, D.; He, X.; Han, J.; and Huang, T. S. 2011. Graphregularized nonnegative matrix factorization for data representation. Pattern Analysis and Machine Intelligence, IEEE
Transactions on 33(8):1548–1560.
Cai, D.; He, X.; and Han, J. 2011. Locally consistent concept factorization for document clustering. Knowledge andData Engineering, IEEE Transactions on 23(6):902–913.
Dhillon, I. S.; Guan, Y.; and Kulis, B. 2004. Kernel k-means:
spectral clustering and normalized cuts. In Proceedings ofthe tenth ACM SIGKDD international conference on Knowledge discovery and data mining, 551–556. ACM.
Dickinson, P. J., and Gijben, L. 2014. On the computational
complexity of membership problems for the completely positive cone and its dual. Computational Optimization and Applications 57(2):403–415.
Ding, C. H.; He, X.; and Simon, H. D. 2005. On the equivalence of nonnegative matrix factorization and spectral clustering. In SDM, volume 5, 606–610.
Ding, C.; Li, T.; and Peng, W. 2008. On the equivalence
between non-negative matrix factorization and probabilistic
latent semantic indexing. Computational Statistics & Data
Analysis 52(8):3913–3927.
Erhan, D.; Bengio, Y.; Courville, A.; Manzagol, P.-A.; Vincent, P.; and Bengio, S. 2010. Why does unsupervised pretraining help deep learning? The Journal of Machine Learning Research 11(3):625–660.
He, Z.; Xie, S.; Zdunek, R.; Zhou, G.; and Cichocki, A.
2011. Symmetric nonnegative matrix factorization: Algorithms and applications to probabilistic clustering. Neural
Networks, IEEE Transactions on 22(12):2117–2131.
Jain, A. K., and Dubes, R. C. 1988. Algorithms for clustering data. Prentice-Hall, Inc.
Jain, A. K. 2010. Data clustering: 50 years beyond k-means.
Pattern Recognition Letters 31(8):651–666.
Kim, H., and Park, H. 2008a. Nonnegative matrix factorization based on alternating nonnegativity constrained leastsquares and active set method. SIAM Journal on Matrix
Analysis and Applications 30(2):713–730.
Kim, J., and Park, H. 2008b. Sparse nonnegative matrix
factorization for clustering.
Kim, J., and Park, H. 2008c. Toward faster nonnegative
matrix factorization: A new algorithm and comparisons. In

Data Mining, 2008. ICDM’08. Eighth IEEE International
Conference on, 353–362. IEEE.
Kleinberg, J. 2003. An impossibility theorem for clustering. Advances in neural information processing systems
463–470.
Kuang, D.; Park, H.; and Ding, C. H. 2012. Symmetric nonnegative matrix factorization for graph clustering. In SDM,
volume 12, 106–117.
Lovász, L., and Plummer, M. D. 1986. Matching theory.
New York.
Minka, T. P. 2000. Old and new matrix algebra usefulfor statistics. See www. stat. cmu. edu/minka/papers/matrix.
html.
Murphy, K. P. 2012. Machine learning: a probabilistic perspective. The MIT Press.
Ng, A. Y.; Jordan, M. I.; Weiss, Y.; et al. 2002. On spectralclustering: Analysis and an algorithm. Advances in neural
information processing systems 2:849–856.
Pei, Y., and Zaı̈ane, O. 2006. A synthetic data generator forclustering and outlier analysis. Department of Computing
science, University of Alberta, edmonton, AB, Canada.
Seung, D., and Lee, L. 2001. Algorithms for non-negative
matrix factorization. Advances in neural information processing systems 13:556–562.
Shi, J., and Malik, J. 2000. Normalized cuts and image
segmentation. Pattern Analysis and Machine Intelligence,
IEEE Transactions on 22(8):888–905.
Steinbach, M.; Karypis, G.; Kumar, V.; et al. 2000. A comparison of document clustering techniques. In KDD workshop on text mining, volume 400, 525–526. Boston.
Thalamuthu, A.; Mukhopadhyay, I.; Zheng, X.; and Tseng,
G. C. 2006. Evaluation and comparison of gene clustering methods in microarray analysis. Bioinformatics
22(19):2405–2412.
Vavasis, S. A. 2009. On the complexity of nonnegative matrix factorization. SIAM Journal on Optimization
20(3):1364–1377.
Vogt, J. E.; Prabhakaran, S.; Fuchs, T. J.; and Roth, V. 2010.
The translation-invariant wishart-dirichlet process for clustering distance data. In Proceedings of the 27th International
Conference on Machine Learning (ICML-10), 1111–1118.
Wilcoxon, F., and Wilcox, R. A. 1964. Some rapid approximate statistical procedures. Lederle Laboratories.
Xu, W.; Liu, X.; and Gong, Y. 2003. Document clustering based on non-negative matrix factorization. In Proceedings of the 26th annual international ACM SIGIR conference
on Research and development in informaion retrieval, 267–
273. ACM.
Zelnik-Manor, L., and Perona, P. 2004. Self-tuning spectralclustering. In Advances in neural information processing
systems, 1601–1608.