I have some points $X=\{x_1,...,x_n\}$ in $R^p$, and I want to cluster the points so that:
Each cluster contains an equal number of elements of $X$. (Assume that the number of clusters divides $n$.)
Each cluster is "spatially cohesive" in some sense, like the clusters from $k$-means.
It's easy to think of a lot of clustering procedures that satisfy one or the other of these, but does anyone know of a way to get both at once?
I suggest a two-step approach:
get a good initial estimates of the cluster centers, e.g. using hard or fuzzy K-means.
Use Global Nearest Neighbor assignment to associate points with cluster centers: Calculate a distance matrix between each point and each cluster center (you can make the problem a bit smaller by only calculating reasonable distances), replicate each cluster center X times, and solve the linear assignment problem. You'll get, for each cluster center, exactly X matches to data points, so that, globally, the distance between data points and cluster centers is minimized.
Note that you can update cluster centers after step 2 and repeat step 2 to basically run K-means with fixed number of points per cluster. Still, it will be a good idea to get a good initial guess first.
Try this k-means variation:
Initialization:
k centers from the dataset at random, or even better using kmeans++ strategyIn the end, you should have a paritioning that satisfies your requirements of the +-1 same number of objects per cluster (make sure the last few clusters also have the right number. The first m clusters should have ceil objects, the remainder exactly floor objects.)
Iteration step:
Requisites: a list for each cluster with "swap proposals" (objects that would prefer to be in a different cluster).
E step: compute the updated cluster centers as in regular k-means
M step: Iterating through all points (either just one, or all in one batch)
Compute nearest cluster center to object / all cluster centers that are closer than the current clusters. If it is a different cluster:
The cluster sizes remain invariant (+- the ceil/floor difference), an objects are only moved from one cluster to another as long as it results in an improvement of the estimation. It should therefore converge at some point like k-means. It might be a bit slower (i.e. more iterations) though.
I do not know if this has been published or implemented before. It's just what I would try (if I would try k-means. there are much better clustering algorithms.)
A good place to start might be with the k-means implementation in ELKI, which already seems to support three different initializations (including k-means++), and the authors said they also want to have different iteration strategys, to cover all the various common variants in a modular fashion (e.g. Lloyd, MacQueen, ...).
Recently I needed this myself for a not very large dataset. My answer, although it has a relatively long running time, is guaranteed to converge to a local optimum.
def eqsc(X, K=None, G=None): "equal-size clustering based on data exchanges between pairs of clusters" from scipy.spatial.distance import pdist, squareform from matplotlib import pyplot as plt from matplotlib import animation as ani from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection def error(K, m, D): """return average distances between data in one cluster, averaged over all clusters""" E = 0 for k in range(K): i = numpy.where(m == k)[0] # indeces of datapoints belonging to class k E += numpy.mean(D[numpy.meshgrid(i,i)]) return E / K numpy.random.seed(0) # repeatability N, n = X.shape if G is None and K is not None: G = N // K # group size elif K is None and G is not None: K = N // G # number of clusters else: raise Exception('must specify either K or G') D = squareform(pdist(X)) # distance matrix m = numpy.random.permutation(N) % K # initial membership E = error(K, m, D) # visualization #FFMpegWriter = ani.writers['ffmpeg'] #writer = FFMpegWriter(fps=15) #fig = plt.figure() #with writer.saving(fig, "ec.mp4", 100): t = 1 while True: E_p = E for a in range(N): # systematically for b in range(a): m[a], m[b] = m[b], m[a] # exchange membership E_t = error(K, m, D) if E_t < E: E = E_t print("{}: {}<->{} E={}".format(t, a, b, E)) #plt.clf() #for i in range(N): #plt.text(X[i,0], X[i,1], m[i]) #writer.grab_frame() else: m[a], m[b] = m[b], m[a] # put them back if E_p == E: break t += 1 fig, ax = plt.subplots() patches = [] for k in range(K): i = numpy.where(m == k)[0] # indeces of datapoints belonging to class k x = X[i] patches.append(Polygon(x[:,:2], True)) # how to draw this clock-wise? u = numpy.mean(x, 0) plt.text(u[0], u[1], k) p = PatchCollection(patches, alpha=0.5) ax.add_collection(p) plt.show() if __name__ == "__main__": N, n = 100, 2 X = numpy.random.rand(N, n) eqsc(X, G=3) This is an optimisation problem. We have an open source java library which solves this problem (clustering where quantity per cluster must be between set ranges). You'd need your total number of points to be maximum of a few thousand though - no more than 5000 or maybe 10000.
The library is here:
https://github.com/PGWelch/territorium/tree/master/territorium.core
The library itself is setup for geographic / GIS type problems - so you will see references to X and Ys, latitudes and longitudes, customers, distance and time, etc. You can just ignore the 'geographic' elements though and use it as a pure clusterer.
You provide a set of initially empty input clusters each with a min and max target quantity. The clusterer will assign points to your input clusters, using a heuristic-based optimisation algorithm (swaps, moves etc). In the optimisation it firstly prioritises keeping each cluster within its min and max quantity range and then secondly minimises the distances between all points in the cluster and the cluster's central point, so a cluster is spatially cohesive.
You give the solver a metric function (i.e. distance function) between points using this interface:
The metric is actually structured to return both a distance and 'time', because its designed for travel-based geographic problems, but for arbitrary clustering problems just set 'time' to be zero and distance to be your actual metric you're using between points.
You'd setup your problem in this class:
Your points would be the 'Customers' and their quantity would be 1. In the customer class ensure you set costPerUnitTime = 0 and costPerUnitDistance=1 assuming you're returning your metric distance in the 'distance' field returned by the TravelMatrix.
See here for an example of running the solver:
I suggest the recent paper Discriminative Clustering by Regularized Information Maximization (and references therein). Specifically, Section 2 talks about class balance and cluster assumption.
