We show that k-means clustering is a matrix factorization problem. Seen from this point of view, k-means clustering can be computed using alternating least squares techniques and we show how the constrained optimization steps involved in this procedure can be solved e ciently using the Frank-Wolfe algorithm.
In this paper, we are concerned with theoretical aspects of machine learning. In particular, we revisit k-means clustering and investigate it from the point of view of data matrix factorization.
The k-means procedure is a popular technique to cluster a data set X of numerical data into subsets C1; : : : ; Ck. The underlying ideas are intuitive and simple and most properties of k-means clustering are text book material [1, 2]. Adding to this material, several authors have recently argued that k-means clustering can be understood as a constrained matrix factorization problem [3{7]. However, reading the related literature, one cannot but notice that most authors consider this fact self explanatory and appeal to intuition.
Our goals with this paper are therefore as follows: i) we provide a rigorous proof for the equivalence of k-means clustering and constrained data matrix factorization. ii) we show that the matrix factorization point of view immediately reveals several properties of k-means clustering which are usually di cult to work out. In particular, we show that k-means clustering is an integer programming problem and therefore NP-hard, that k-means clustering allows for invoking the kernel trick, and that k-means clustering is closely related to archetypal analysis [8, 9] and non-negative matrix factorization [10, 11]. iii) we show that the matrix factorization perspective leads to yet another algorithm for computing k-means clustering and we discuss how to e ciently implement it using the Frank-Wolfe optimization scheme [12{14].
We begin our presentation with a brief summery of the traditional view on k-means clustering and then move on to the matrix factorization perspective. Our discussion assumes that readers are familiar with theory and practice of matrix factorization for data mining and machine learning. Those interested in a gentle introduction into the underlying mathematical ideas are referred to [11].
k-Means Clustering: Known Properties and Algorithms Given a set X = fx1; : : : ; xng of data points xj 2 Rm, hard k-means clustering attempts to partition the data into k clusters C1; : : : ; Ck such that Ci X, Ci \ Cl = ;, and C1 [ C2 [ : : : [ Ck = X. In particular, hard k-means clustering is a prototype-based clustering technique because it understands clusters to be de ned in terms of prototypes or cluster centroids 1; : : : ; k 2 Rm, namely Ci = nxj 2 X i 2 xj l 2 8 l 6= i o
The problem at the heart of hard k-means clustering is therefore to search for k appropriate cluster centroids which are typically determined as the minimizers of the following objective function
k E(k) = X X i=1 xj2Ci i
k 2 = X n
X zij xj i=1 j=1 i 2 where the zij = are binary indicator variables which indicate whether or not data point xj belongs to cluster Ci.
Since hard k-means clustering aims at disjoint clusters where each xj is assigned to one and only one Ci, we point out the following important properties of the zij 2 f0; 1g. If we x the data index j and sum over the cluster index i, we obtain the number of clusters data point xj is assigned to, namely xj xj k X zij = 1: i=1 n X zij = jCij = ni j=1 Also, by xing the cluster index i and summing over the data index j, we nd where ni indicates the number of data points assigned to cluster Ci.
Although the objective in (
k-Means Clustering Is Data Matrix Factorization
Having recalled properties of hard k-means clustering in the previous section, our
goal is now to rigorously establish that k-means clustering is matrix factorization.
In other words, given the objective in (
X
M Z 2 where
X 2 Rm n is a column matrix of n data vectors xj 2 Rm M 2 Rm k is a column matrix of k cluster centroids i 2 Rm Z 2 Rk n is a matrix of binary indicator variables such that zij = Throughout, we write xj to denote j-th column vector of a matrix X. To refer to the (l; j) element of a matrix X, we either write xlj or X lj . Moreover, subscripts or summation indices i will be understood to range from 1 to k (the number of clusters), subscripts or summation indices j will range from 1 up to n (the number of data), and subscripts or summation indices l will be used to expand inner products between vectors or rows and columns of matrices.
Finally, regarding the squared Frobenius norm of a matrix, we recall that
X 2 = X xl2j = X xj 2 = X xjT xj = X XT X jj = tr XT X (
Step by Step Derivation of (
Z ij ZT ji0 = which is then to say that the matrix ZZT is a diagonal matrix where X (ni; if i = i0 0;
otherwise.
X zij zi0j =
j
Having discussed this property of Z, we are now positioned to establish the
equality in (
X zij xj i;j i 2 = X zij xjT xj 2xjT i + i i
T i;j {z
T1 = X zij xjT xj 2 X zij xjT i + X zij iT i : i;j i;j i;j | } | } | } {z T2
This expansion leads to further insights, if we examine the three terms T1, T2, and T3 one by one. First of all, we nd
T1 = X zij xjT xj = X zij xj 2 = X
i;j i;j j
i
xj 2 X zij = X
xj
2
= tr XT X
where we made use of (
T2 = X zij xj i = X zij
T i;j i;j {z T3 j X xlj li
l = X xlj X li zij
j;l i = X xlj M Z lj
j;l = X
XT jl
M Z lj = X
XT M Z jj = tr XT M Z j;l j i i 2 = X i 2 X zij = X i 2
ni
j
i
(
T3 = X zij i i = X zij
T
i;j i;j
where we applied (
Step 2: Expanding the expression on the right of (
M Z 2 = trh X
M Z T X
M Z i = tr XT X | {Tz4 } 2 tr XT M Z + tr ZT M T M Z | {Tz5 } | {Tz6 } (23)
Given our earlier results, we immediately recognize that T1 = T4 and T2 = T5.
Thus, to establish that (
Regarding T6, we note that, because of the cyclic permutation invariance of the trace operator, we have
tr ZT M T M Z
= tr M T M ZZT : tr M T M ZZT
M T M ZZT
ii = X
i = X X
i = X
i = X i l M T M i 2
ni M T M
il ZZT li
ii ZZT ii
where we used the fact that ZZT is diagonal. This result, however, shows that
T3 = T6 and, consequently, that (
Having proved our central claim in (
k-Means Clustering Is NP Hard
Given the above result, further insights into the nature of k-means clustering
result from eliminating matrix M from the right hand side of (
X
M Z 2 = = 2 M ZZT
XZT which, upon equating to 0, leads to (31) (32) which beautifully re ects the fact that each of the k-means cluster centroids coincides with the mean of the corresponding cluster Ci, namely i
Given the result in (30), we therefore nd that the k-means objective can be cast solely in terms of the data matrix X and the indicator matrix Z i =
Pj zij xj
Pj zij =
Looking at (32), we recognize k-means clustering as an integer programming problem, since it corresponds to the discrete optimization problem of nding a column stochastic binary matrix Z that minimizes the objective in (32). Integer programming problems are NP hard and we can actually read this o (32). Z is an k n binary matrix such that each column contains a single 1; thus, for each column there are k ways to place that 1 and since there are n columns, there are O kn matrices among which we have to determine the optimal one. Observe that the squared Frobenius norm in (32) can also be expanded in terms of trace operators. If we substitute = ZT ZZT 1Z for brevity, we nd X
X 2 = trhXT Xi 2 trhXT X i + trh
XT X i (33) and recognize that each occurrence of the data vectors xi is in form of an inner product, because XT X ij = xiT xj .
This, however, shows that k-means allows for invoking the kernel trick as we may replace the n n Gramian XT X by an n n kernel matrix K whose entries correspond to kernel evaluations k(xi; xj ). 4.3 k-Means Clustering, Archetypal Analysis and Non-Negative Matrix Factorization Even further insights into the nature of k-means clustering arise, if we substitute Y = ZT ZZT 1 so that we can write M = XY . This again reveals that the centroid vectors i, i.e. the columns of M , are convex combinations of data points xj . That is, i = Xyi where yi is an n dimensional vector with ni entries equal to 1=ni and n ni entries equal to 0.
We note that, by de nition, Y is a column stochastic matrix. It is nonnegative, i.e. Y 0, and its columns sum to one, i.e. Pj yji = 1T yi = 1. Moreover, each column will have high entropy H(yi) = Pj yji log yji 0.
We also recall that, as a binary matrix, matrix Z will obey Z 0, 1T zj = 1, and H(zj ) = Pi zij log zij = 0.
Hard k-Means Clustering Given these prerequisites, we can therefore express hard k-means clustering as a constrained quadratic optimization problem Archetypal Analysis If we drop the entropy constraints in (34), we obtain min Y ;Z and recover a problem known as archetypal analysis (AA) [8, 9]. Dropping the entropy constraints has an interesting e ect. Instead of computing basis vectors M = XY that correspond to local means, AA determines basis vectors that are extreme points of the data in X. In fact, the archetypes in matrix M reside on the data convex hull [22].
Non-Negative Matrix Factorization If we further drop the sum-to-one constraints in (35), we obtain min Y ;Z s.t.
X
XY Z
2 Y Z 0 0 a problem known as non-negative matrix factorization (NMF) [10]. Dropping the stochastic constraints has the e ect that NMF computes basis vectors M = XY that are conic combinations of the data in X.
The expressions in (34), (35), and (36) therefore reveal k-means clustering to be a particularly severely constrained quadratic optimization problem. This is interesting in so far as algorithms for computing k-means clustering conceptually much simpler than AA or NMF algorithms. Nevertheless, it now appears as if methods that have been developed for AA and NMF might also apply to k-means clustering. In the next section, we show that this is indeed the case.
Since we just found that hard k-means clustering is indeed a constrained form of archetypal analysis, the question is if algorithms that have been developed for archetypal analysis can be used to compute k-means clustering.
In order to see how this can be accomplished, we note that the objective function in (34) is convex in either Y or Z but not in the their product Y Z. An idea for solving the problem in (34) could therefore be to apply the following constrained alternating least squares procedure 1) randomly initialize Y and Z under the appropriate constraints 2) x matrix Z and update Y to the solution of 3) x matrix Y and update Z to the solution of
XY Z XY Z 2 2 min f (x) x
Second of all, the column entropy H(yi) = Pj yji log yji is a concave function so that H(yi) is convex and we are interested in solutions for Y that maximize H(yi). We can thus rewrite the rst problem in the above procedure as 4) if not converged, continue at 2) for which we note that the seemingly di cult constrained quadratic optimization problems it involves can actually be solved rather easily.
First of all, we observe that the stochastic constraints yi 0 and 1T yi = 1 require the columns of matrix Y to reside in the standard n-simplex (37) (38) (40) min Y s.t.
X
XY Z yi 2 n 1: 2
X H(yi) i
Written in this form, we recognize the problem as a convex minimization problem over a compact convex set. This, however, is to say that it can be tackled using the e cient Frank-Wolfe procedure [12]. The Frank-Wolfe algorithm shown in Alg. 1 solves problems of the form where S Rm is a compact, convex set and f : S ! R is a convex function. The key idea is to compute s 2 S that minimizes sT rf (xt) and to use sub-gradient updates s xt which guarantee that the updates will never leave the feasible set. The e ciency of the algorithm stems from the fact that it turns a quadratic optimization problem into a series of linear optimization problems and we point out that the minimum of a linear function sT rf (x) over a compact convex set will be attained at a vertex of that set.
With respect to our problem of a matrix minimization problem, we note that the gradient rf we are concerned with is given by rY
X
XY Z 2
X H(yi) i
! where the components of matrix L amount to = 2 hXT XY ZZT
XT XZT i
L (41) (39) (42) Algorithm 1 Frank-Wolfe algorithm to solve problem such as in (40) guess a feasible point x0 for t = 0; : : : ; tmax do determine st by solving min sT rf (xt) s2S 2 update the learning rate t = t+2 update the current estimate xt+1 = xt + t (st xt) Lji =
X yji log yji = log yji + 1 @yji j We can use the columns of this gradient matrix G to update the columns yi of Y according to the Frank-Wolfe algorithm. We also point out, the the compact convex set we are dealing with is the standard simplex n 1 whose vertices are given by the standard basis vectors el 2 Rn. The problem of nding the minimizer s 2 n 1 thus simpli es to nding the basis vector el that minimizes elT gi which is simply to determine the minimal entry gli in each column gi of G. All in all, these consideration then lead to Alg. 2 for computing updates of matrix Y .
Similar considerations apply to the problem of computing the updates of the indicator matrix Z and correspondingly lead to Alg. 3.
To conclude this discussion, we point out that the Frank-Wolfe procedure quickly achieves -approximations of the optimal solution that are provably sparse [13]. In fact, one can show that after t iterations the current estimate is O(1=t) from the optimal solution [13] which provides a convenient criterion for choosing the number tmax of iterations to be performed. For further details on the Frank-Wolfe algorithms as well as for a recent excellent survey of projection-free convex optimization over compact convex sets, we refer to [14]. Algorithm 2 Frank-Wolfe procedure to compute Y whose columns are in Require: data matrix X, indicator matrix Z, and parameter tmax 2 N Y [e1; e1; : : : ; e1] where e1 = [1; 0; : : : ; 0]T 2 Rn // initialize n t 0 repeat
G = 2 XT XY ZZT XT XZT L for i = 1; : : : ; k do i0 = argminl Gli zi zi + 2=(t + 2) (ei0 zi) t t + 1 until updates \become small" or t = tmax // compute gradient matrix // update columns yi of Y Algorithm 3 Frank-Wolfe procedure to compute Z whose columns are in Require: data matrix X, coe cient matrix Y , and parameter tmax 2 N Z [e1; e1; : : : ; e1] where e1 = [1; 0; : : : ; 0]T 2 Rk // initialize k t 0 repeat
G = 2 Y T XT XY Z Y T XT X for j = 1; : : : ; n do j0 = argminl Glj zj zj + 2=(t + 2) (ej0 zj ) t t + 1 until updates \become small" or t = tmax // compute gradient matrix // update columns zj of Z k matrix Y n matrix Z 6
In this paper, we were concerned with machine learning theory. In particular, we were concerned with theoretical aspects of k-means clustering.
First of all, we rigorously established that k-means clustering is a constrained data matrix factorization problem. Second of all, this insight allowed us to easily uncover several properties of k-means clustering that are otherwise more di cult to show. Third of all, the matrix factorization point of view on k-means clustering allowed us to reveal its connections to archetypal analysis and non-negative matrix factorization. Finally, given that k-means clustering can be understood as a constrained version of archetypal analysis, we discussed yet another algorithm for k-means clustering. Archetypal analysis is often computed using alternating least squares optimization and we showed how to adapt this idea to k-means clustering. In particular, we discussed that the seemingly di cult constrained optimization problems involved in this procedure can be solved using the e cient Frank-Wolfe procedure for convex optimization over compact convex sets.
Again, the work reported here is of mainly theoretical interest. What is particularly striking is that it established k-means clustering as a more constrained and thus more di cult problem than archetypal analysis or non-negative matrix factorization. Yet, at the same time, traditional algorithms for k-means clustering are considerably simpler than those for the latter problems. This can be seen as a call to arms for it suggests that there may be simpler algorithms for these kind of problems as well. Indeed, techniques such as k-maxoids clustering [23] which were derived from k-means clustering indicate that, say, archetypal analysis should be solvable by simple algorithms, too. 21. Bottou, L., Bengio, Y.: Convergence Properties of the K-Means Algorithms. In:
Proc. NIPS. (1995) 22. Bauckhage, C.: A Note on Archetypal Analysis and the Approximation of Convex
Hulls. arXiv:1410.0642 [cs.NA] (2014) 23. Bauckhage, C., Sifa, R.: k-Maxoids Clustering. In: Proc. KDML-LWA. (2015)