k-means clustering algorithm implementation 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 &lt; 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?
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 &lt; 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")