k-means clustering algorithm implementation

Question

Here is my personal implementation of the clustering k-means algorithm.

from scipy.spatial import distance
import numpy as np
import random
# (x,y) coordinates of a point
X = 0
Y = 1
def get_first(k, points):
 return points[0:k]
def cost(cetroids, clusters):
 cost = 0
 for i in range(len(centroids)):
 centroid = centroids[i]
 cluster = clusters[i]
 for point in cluster:
 cost += (distance.euclidean(centroid, point))**2
 return cost
def compute_centroids(clusters):
 centroids = []
 for cluster in clusters:
 centroids.append(np.mean(cluster, axis=0))
 return centroids
def kmeans(k, centroids, points, method, iter):
 clusters = [[] for i in range(k)]
 for i in range(len(points)):
 point = points[i]
 belongs_to_cluster = closest_centroid(point, centroids)
 clusters[belongs_to_cluster].append(point)
 new_centroids = compute_centroids(clusters)
 if not equals(centroids, new_centroids):
 print "Iteration " + str(iter) + ". Cost [k=" + str(k) + ", " + method + "] = " + str(cost(new_centroids, clusters))
 clusters = kmeans(k, new_centroids, points, method, iter+1)
 return clusters
def closest_centroid(point, centroids):
 min_distance = float('inf')
 belongs_to_cluster = None
 for j in range(len(centroids)):
 centroid = centroids[j]
 dist = distance.euclidean(point, centroid)
 if dist < min_distance:
 min_distance = dist
 belongs_to_cluster = j
 return belongs_to_cluster
def contains(point1, points):
 for point2 in points:
 if point1[X] == point2[X] and point1[Y] == point2[Y]:
 return True
 return False
def equals(points1, points2):
 if len(points1) != len(points2):
 return False
 for i in range (len(points1)):
 point1 = points1[i]
 point2 = points2[i]
 if point1[X] != point2[X] or point1[Y] != point2[Y]:
 return False
 return True
if __name__ == "__main__":
 data = [[-19.0748, -8.536 ],
 [ 22.0108, -10.9737 ],
 [ 12.6597, 19.2601 ],
 [ 11.26884087, 19.90132146 ],
 [ 15.44640731, 21.13121676 ],
 [-20.03865146, -8.820872829],
 [-19.65417726, -8.211477352],
 [-15.97295894, -9.648002534],
 [-18.74359696, -5.383551586],
 [-19.453215, -8.146120006],
 [-16.43074088, -7.524968005],
 [-19.75512437, -8.533215751],
 [-19.56237082, -8.798668569],
 [-19.47135573, -8.057217004],
 [-18.60946986, -4.475888949],
 [-21.59368337, -10.38712463],
 [-15.39158057, -3.8336522 ],
 [-40.0, 40.0 ]]
 k = 4
 # k-means picking the first k points as centroids
 centroids = get_first(k, data)
 clusters = kmeans(k, centroids, data, "first", 1)

I understood and followed the theory of the algorithm, but as you can see, when running the code, the cost on each iteration of the algorithm increases. I am new to python and I think the problem relies in some misunderstanding from my side regarding list manipulation. Any ideas?

• Phrancis: Welcome to Code Review! I hope you get some good answers.
• Daniyal Shahrokhian: Also, I am open to suggestions regarding the quality of my code.

Answer 1

N.b. I'm pretty sure you can find optimised replacements of many of these functions in either the NumPy or SciPy library.

Edit: I don't immediately see a reason for a slowdown except the addition of new clusters. A bigger sample set that shows the behaviour would be nice, together with some timings and probably running a profiler on it.

In general you should likely use a NumPy array or matrix instead of lists of lists as they are optimised for storing numeric data(!). That would likely also eliminate the need for some of these functions or considerably reduce their length. Also store all data in the same containers - that way you don't have to create new functions like equals and contains yourself to compare between different representations.

• The X/Y definitions at the start are more of a WTF for me. I'd get rid of that by writing code that's either independent on the number of dimensions, or just use 0/1 - IMO it's not magic enough to warrant separate constants.
• get_first isn't neccessary - the pattern is obvious enough to not require encapsulation in a function. The zero is also optional.
• Typo in cost signature, should be centroids.
• The pattern for i in range(len(...)): appears a couple of times and should be replaced with proper iteration over some helper generator. E.g. in cost the parallel iteration over both centroids and clusters should be done by zip (itertools.izip in Python 2.7 if used in the same way). It can also be simplified by using another SciPy function, cdist.
• If you still need the index, use enumerate.
• You can also often get away with using the squared euclidean if the exact value isn't relevant, just the relation between different distances.
• The pattern to initialise a list of empty lists should probably use _ instead of i to make it obvious that the variable serves to further purpose.
• print should be used as a function to make it more uniform. Also, use one of the various formatting options (% or format) to get rid of the additional str calls.
• compute_centroids can be simplified with a list comprehension.

Looks like this now:

from scipy.spatial import distance
import numpy as np
import random
from itertools import izip
def cost(centroids, clusters):
 return sum(distance.cdist([centroid], cluster, 'sqeuclidean').sum()
 for centroid, cluster in izip(centroids, clusters))
def compute_centroids(clusters):
 return [np.mean(cluster, axis=0) for cluster in clusters]
def kmeans(k, centroids, points, method):
 clusters = [[] for _ in range(k)]
 for point in points:
 clusters[closest_centroid(point, centroids)].append(point)
 new_centroids = compute_centroids(clusters)
 if not equals(centroids, new_centroids):
 print("cost [k={}, {}] = {}".format(k, method, cost(new_centroids, clusters)))
 clusters = kmeans(k, new_centroids, points, method)
 return clusters
def closest_centroid(point, centroids):
 min_distance = float('inf')
 belongs_to_cluster = None
 for j, centroid in enumerate(centroids):
 dist = distance.sqeuclidean(point, centroid)
 if dist < min_distance:
 min_distance = dist
 belongs_to_cluster = j
 return belongs_to_cluster
def contains(point1, points):
 for point2 in points:
 if point1[0] == point2[0] and point1[1] == point2[1]:
 # if all(x == y for x, y in izip(points1, points2)):
 return True
 return False
def equals(points1, points2):
 if len(points1) != len(points2):
 return False
 for point1, point2 in izip(points1, points2):
 if point1[0] != point2[0] or point1[1] != point2[1]:
 # if any(x != y for x, y in izip(points1, points2)):
 return False
 return True
if __name__ == "__main__":
 data = [[-19.0748, -8.536 ],
 [ 22.0108, -10.9737 ],
 [ 12.6597, 19.2601 ],
 [ 11.26884087, 19.90132146 ],
 [ 15.44640731, 21.13121676 ],
 [-20.03865146, -8.820872829],
 [-19.65417726, -8.211477352],
 [-15.97295894, -9.648002534],
 [-18.74359696, -5.383551586],
 [-19.453215, -8.146120006],
 [-16.43074088, -7.524968005],
 [-19.75512437, -8.533215751],
 [-19.56237082, -8.798668569],
 [-19.47135573, -8.057217004],
 [-18.60946986, -4.475888949],
 [-21.59368337, -10.38712463],
 [-15.39158057, -3.8336522 ],
 [-40.0, 40.0 ]]
 k = 4
 # k-means picking the first k points as centroids
 centroids = data[:k]
 clusters = kmeans(k, centroids, data, "first")

• Daniyal Shahrokhian: This is... just prodigious. Sharp, neat and complete answer. I will revise the whole thing tomorrow morning. Thanks a lot @ferada.