Voici un petit kmeans qui utilise n'importe quelle distance parmi la vingtaine de distances dans scipy.spatial.distance ou une fonction utilisateur.
Les commentaires sont les bienvenus (il n'y a eu qu'un seul utilisateur jusqu'à présent, ce qui n'est pas suffisant) ; En particulier, quels sont vos N, dim, k, métriques ?
#!/usr/bin/env python
# kmeans.py using any of the 20-odd metrics in scipy.spatial.distance
# kmeanssample 2 pass, first sample sqrt(N)
from __future__ import division
import random
import numpy as np
from scipy.spatial.distance import cdist # $scipy/spatial/distance.py
# http://docs.scipy.org/doc/scipy/reference/spatial.html
from scipy.sparse import issparse # $scipy/sparse/csr.py
__date__ = "2011-11-17 Nov denis"
# X sparse, any cdist metric: real app ?
# centres get dense rapidly, metrics in high dim hit distance whiteout
# vs unsupervised / semi-supervised svm
#...............................................................................
def kmeans( X, centres, delta=.001, maxiter=10, metric="euclidean", p=2, verbose=1 ):
""" centres, Xtocentre, distances = kmeans( X, initial centres ... )
in:
X N x dim may be sparse
centres k x dim: initial centres, e.g. random.sample( X, k )
delta: relative error, iterate until the average distance to centres
is within delta of the previous average distance
maxiter
metric: any of the 20-odd in scipy.spatial.distance
"chebyshev" = max, "cityblock" = L1, "minkowski" with p=
or a function( Xvec, centrevec ), e.g. Lqmetric below
p: for minkowski metric -- local mod cdist for 0 < p < 1 too
verbose: 0 silent, 2 prints running distances
out:
centres, k x dim
Xtocentre: each X -> its nearest centre, ints N -> k
distances, N
see also: kmeanssample below, class Kmeans below.
"""
if not issparse(X):
X = np.asanyarray(X) # ?
centres = centres.todense() if issparse(centres) \
else centres.copy()
N, dim = X.shape
k, cdim = centres.shape
if dim != cdim:
raise ValueError( "kmeans: X %s and centres %s must have the same number of columns" % (
X.shape, centres.shape ))
if verbose:
print "kmeans: X %s centres %s delta=%.2g maxiter=%d metric=%s" % (
X.shape, centres.shape, delta, maxiter, metric)
allx = np.arange(N)
prevdist = 0
for jiter in range( 1, maxiter+1 ):
D = cdist_sparse( X, centres, metric=metric, p=p ) # |X| x |centres|
xtoc = D.argmin(axis=1) # X -> nearest centre
distances = D[allx,xtoc]
avdist = distances.mean() # median ?
if verbose >= 2:
print "kmeans: av |X - nearest centre| = %.4g" % avdist
if (1 - delta) * prevdist <= avdist <= prevdist \
or jiter == maxiter:
break
prevdist = avdist
for jc in range(k): # (1 pass in C)
c = np.where( xtoc == jc )[0]
if len(c) > 0:
centres[jc] = X[c].mean( axis=0 )
if verbose:
print "kmeans: %d iterations cluster sizes:" % jiter, np.bincount(xtoc)
if verbose >= 2:
r50 = np.zeros(k)
r90 = np.zeros(k)
for j in range(k):
dist = distances[ xtoc == j ]
if len(dist) > 0:
r50[j], r90[j] = np.percentile( dist, (50, 90) )
print "kmeans: cluster 50 % radius", r50.astype(int)
print "kmeans: cluster 90 % radius", r90.astype(int)
# scale L1 / dim, L2 / sqrt(dim) ?
return centres, xtoc, distances
#...............................................................................
def kmeanssample( X, k, nsample=0, **kwargs ):
""" 2-pass kmeans, fast for large N:
1) kmeans a random sample of nsample ~ sqrt(N) from X
2) full kmeans, starting from those centres
"""
# merge w kmeans ? mttiw
# v large N: sample N^1/2, N^1/2 of that
# seed like sklearn ?
N, dim = X.shape
if nsample == 0:
nsample = max( 2*np.sqrt(N), 10*k )
Xsample = randomsample( X, int(nsample) )
pass1centres = randomsample( X, int(k) )
samplecentres = kmeans( Xsample, pass1centres, **kwargs )[0]
return kmeans( X, samplecentres, **kwargs )
def cdist_sparse( X, Y, **kwargs ):
""" -> |X| x |Y| cdist array, any cdist metric
X or Y may be sparse -- best csr
"""
# todense row at a time, v slow if both v sparse
sxy = 2*issparse(X) + issparse(Y)
if sxy == 0:
return cdist( X, Y, **kwargs )
d = np.empty( (X.shape[0], Y.shape[0]), np.float64 )
if sxy == 2:
for j, x in enumerate(X):
d[j] = cdist( x.todense(), Y, **kwargs ) [0]
elif sxy == 1:
for k, y in enumerate(Y):
d[:,k] = cdist( X, y.todense(), **kwargs ) [0]
else:
for j, x in enumerate(X):
for k, y in enumerate(Y):
d[j,k] = cdist( x.todense(), y.todense(), **kwargs ) [0]
return d
def randomsample( X, n ):
""" random.sample of the rows of X
X may be sparse -- best csr
"""
sampleix = random.sample( xrange( X.shape[0] ), int(n) )
return X[sampleix]
def nearestcentres( X, centres, metric="euclidean", p=2 ):
""" each X -> nearest centre, any metric
euclidean2 (~ withinss) is more sensitive to outliers,
cityblock (manhattan, L1) less sensitive
"""
D = cdist( X, centres, metric=metric, p=p ) # |X| x |centres|
return D.argmin(axis=1)
def Lqmetric( x, y=None, q=.5 ):
# yes a metric, may increase weight of near matches; see ...
return (np.abs(x - y) ** q) .mean() if y is not None \
else (np.abs(x) ** q) .mean()
#...............................................................................
class Kmeans:
""" km = Kmeans( X, k= or centres=, ... )
in: either initial centres= for kmeans
or k= [nsample=] for kmeanssample
out: km.centres, km.Xtocentre, km.distances
iterator:
for jcentre, J in km:
clustercentre = centres[jcentre]
J indexes e.g. X[J], classes[J]
"""
def __init__( self, X, k=0, centres=None, nsample=0, **kwargs ):
self.X = X
if centres is None:
self.centres, self.Xtocentre, self.distances = kmeanssample(
X, k=k, nsample=nsample, **kwargs )
else:
self.centres, self.Xtocentre, self.distances = kmeans(
X, centres, **kwargs )
def __iter__(self):
for jc in range(len(self.centres)):
yield jc, (self.Xtocentre == jc)
#...............................................................................
if __name__ == "__main__":
import random
import sys
from time import time
N = 10000
dim = 10
ncluster = 10
kmsample = 100 # 0: random centres, > 0: kmeanssample
kmdelta = .001
kmiter = 10
metric = "cityblock" # "chebyshev" = max, "cityblock" L1, Lqmetric
seed = 1
exec( "\n".join( sys.argv[1:] )) # run this.py N= ...
np.set_printoptions( 1, threshold=200, edgeitems=5, suppress=True )
np.random.seed(seed)
random.seed(seed)
print "N %d dim %d ncluster %d kmsample %d metric %s" % (
N, dim, ncluster, kmsample, metric)
X = np.random.exponential( size=(N,dim) )
# cf scikits-learn datasets/
t0 = time()
if kmsample > 0:
centres, xtoc, dist = kmeanssample( X, ncluster, nsample=kmsample,
delta=kmdelta, maxiter=kmiter, metric=metric, verbose=2 )
else:
randomcentres = randomsample( X, ncluster )
centres, xtoc, dist = kmeans( X, randomcentres,
delta=kmdelta, maxiter=kmiter, metric=metric, verbose=2 )
print "%.0f msec" % ((time() - t0) * 1000)
# also ~/py/np/kmeans/test-kmeans.py
Quelques notes ajoutées le 26 mars 2012 :
1) pour la distance cosinus, il faut d'abord normaliser tous les vecteurs de données à |X| = 1 ; puis
cosinedistance( X, Y ) = 1 - X . Y = Euclidean distance |X - Y|^2 / 2
est rapide. Pour les vecteurs de bits, gardez les normes séparément des vecteurs au lieu de les étendre aux flottants (bien que certains programmes puissent faire l'expansion pour vous). Pour les vecteurs épars, disons 1 % de N, X . Y devrait prendre un temps O( 2 % N ), espace O(N) ; mais je ne sais pas quels programmes font cela.
2) Scikit-learn clustering donne une excellente vue d'ensemble des k-means, mini-batch-k-means ... avec du code qui fonctionne sur les matrices scipy.sparse.
3) Vérifiez toujours la taille des clusters après le k-means. Si vous vous attendez à des clusters de taille à peu près égale, mais qu'ils sont en réalité [44 37 9 5 5] %
... (bruit de grincement de tête).
44 votes
Notez que k-means est conçu pour la distance euclidienne. . Elle peut cesser de converger avec d'autres distances, lorsque le moyenne n'est plus la meilleure estimation du "centre" de la grappe.
2 votes
Pourquoi k-means ne fonctionne que pour la distance euclidienne ?
9 votes
@Anony-Mousse Il est incorrect de dire que k-means est uniquement conçu pour la distance euclidienne. Il peut être modifié pour fonctionner avec toute métrique de distance valide définie sur l'espace d'observation. Par exemple, regardez l'article sur k-medoids .
2 votes
PAM (alias k-medoids) est un algorithme très différent. Il est apparenté aux k-means mais beaucoup plus chers.
5 votes
@curious : le moyenne minimise les différences au carré (= distance euclidienne au carré). Si vous voulez une fonction de distance différente, vous devez remplacer l'option moyenne avec une estimation appropriée du centre. Les K-médoïdes sont un tel algorithme, mais trouver le médoïde est beaucoup plus coûteux.
4 votes
Quelque peu pertinent ici : il y a actuellement un demande de retrait ouverte mettant en œuvre le noyau K-Means. Lorsqu'il sera terminé, vous pourrez spécifier votre propre noyau pour le calcul.
1 votes
@ely. "Il est incorrect de dire que k-means est uniquement conçu pour la distance euclidienne." Non, ce n'est pas incorrect, IMHO. K-means et K-medoids peuvent être liés, mais ce sont des algorithmes différents avec des modèles mathématiques sous-jacents différents, et donc des conditions de convergence différentes. K-means suppose une distance euclidienne. K-medoids suppose une distance de Manhattan. Veuillez me corriger si je me trompe.
0 votes
@ChirazBenAbdelkader Il s'agit du même algorithme avec spécifiquement le même modèle sous-jacent. Ils ne diffèrent que par le calcul spécifique de l'exemplaire utilisé (qu'il s'agisse d'un centroïde de groupe ou d'un medoid de groupe réel). K-means fait référence à un famille d'algorithmes qui utilisent tous le même modèle sous-jacent, mais avec des notions différentes de distance ou d'exemplarité.
0 votes
@ely. Je suis partiellement d'accord avec vous. Peut-être que je coupe les cheveux en quatre. Mais cela dépend vraiment de ce que vous considérez comme le " même " modèle. Oui, Kmeans et Kmedoids sont basés sur le même modèle générique. Mais ils sont suffisamment différents et ne sont certainement PAS interchangeables dans la pratique.
1 votes
@ChirazBenAbdelkader De nombreux algorithmes généraux s'accompagnent de variations spécifiques. Par exemple, l'"algorithme" de SVM serait strictement différent, dans un sens pédant, si vous utilisez un noyau RBF ou un noyau polynomial, etc., et les deux choses ne seraient certainement pas facilement interchangeables dans la pratique. Mais il serait stupide de dire que le SVM avec un noyau RBF est "complètement différent" de celui avec un noyau polynomial. Il s'agit clairement du même algorithme, mais un sous-ensemble de l'algorithme peut être interchangé comme hyperparamètre. C'est la même chose avec les algorithmes k-means.
1 votes
Par exemple, envisagez d'utiliser simplement différents noyaux de dissimilarité pour un ensemble donné de points de données, comme dans l'outil de scipy
pdist
. Vous pouvez également envisager de minimiser la norme L1 ou la divergence KL si vous disposez de points de données soumis à une contrainte d'éparpillement ou qui sont des distributions de probabilité. Il n'y a aucune raison de ne pas exécuter l'algorithme k-means sur ces types de données, et d'utiliser simplement une "distance entre les points" différente et appropriée pour minimiser la fonction de perte par rapport aux centres candidats.