k-means#

Dénomination française : algorithme des centres mobiles.

Principe#

Les centres mobiles ou nuées dynamiques sont un algorithme de classification non supervisée. A partir d’un ensemble de points, il détermine pour un nombre de classes fixé, une répartition des points qui minimise un critère appelé inertie ou variance intra-classe.

Algorithme A1 : centre mobile, k-means

On considère un ensemble de points :

\left(X_i\right)_{1\leqslant i\leqslant P}\in\left(\R^N\right)^P

A chaque point est associée une classe : \left(c_i\right)_{1\leqslant i\leqslant P}\in\left\{1,...,C\right\}^P. On définit les barycentres des classes : \left( G_i\right)_{1\leqslant i\leqslant C}\in\left(\R^N\right)^C.

Initialisation

L’initialisation consiste à choisir pour chaque point une classe aléatoirement dans \left\{1,...,C\right\}. On pose t = 0.

Calcul des barycentres

for k in 1..C
G_k^t \longleftarrow \sum_{i=1}^P X_i \, \mathbf{1}_{\left\{c_i^t=k\right\}} \sum_{i=1}^P \mathbf{1}_{\left\{c_i^t=k\right\}}

Calcul de l’inertie

\begin{array}{lll} I^t &\longleftarrow& \sum_{i=1}^P \; d^2\left(X_i, G_{c_i^t}^t\right) \\ t &\longleftarrow& t+1 \end{array}

if t > 0 et I_t \sim I_{t-1}
arrêt de l’algorithme

Attribution des classes

for in 1..P
c_i^{t+1} \longleftarrow \underset{k}{\arg\min} \; d\left( X_{i},G_{k}^{t}\right)
où d\left(X_i,G_k^t\right) est la distance entre X_i et G_k^t

Retour à l’étape du calcul des barycentres jusqu’à convergence de l’inertie I^t.

Théorème T1 : convergence des k-means

Quelque soit l’initialisation choisie, la suite \pa{I_t}_{t\supegal 0} construite par l’algorithme des k-means converge.

La démonstration du théorème nécessite le lemme suivant.

Lemme L1 : inertie minimum

Soit \vecteur{X_1}{X_P} \in \pa{\R^N}^P, P points de \R^N, le minimum de la quantité Q\pa{Y \in \R^N} :

\begin{eqnarray} Q\pa{Y} &=& \sum_{i=1}^P \; d^2\pa{X_i,Y} \end{eqnarray}

est atteint pour Y=G=\dfrac{1}{P} \sum_{i=1}^{P} X_i le barycentre des points \vecteur{X_1}{X_P}.

Soit \vecteur{X_1}{X_P} \in \pa{\R^N}^P, P points de \R^N.

\begin{eqnarray*} \sum_{i=1}^{P} \overrightarrow{GX_{i}} = \overrightarrow{0} &\Longrightarrow& \sum_{i=1}^{P} d^2\pa{X_i,Y} = \sum_{i=1}^{P} d^2\pa{X_i,G}+ P \, d^2\pa{G,Y} \\ &\Longrightarrow& \underset{Y\in\R^N}{\arg\min} \; \sum_{i=1}^{P} d^2\pa{X_i,Y} = \acc{G} \end{eqnarray*}

On peut maintenant démontrer le théorème. L’étape d’attribution des classes consiste à attribuer à chaque point le barycentre le plus proche. On définit J_t par :

\begin{eqnarray} J^{t+1} &=& \sum_{i=1}^{P} \; d^2\pa{ X_i, G_{c_i^{t+1}}^t} \end{eqnarray}

On en déduit que :

\begin{eqnarray} J^{t+1} &=& \sum_{i, c_i^t \neq c_i^{t+1}} \; d^2\pa{ X_i, G_{c_i^{t+1}}^t} + J^{t+1} \sum_{i, c_i^t = c_i^{t+1}} \; d^2\pa{ X_i, G_{c_i^{t+1}}^t} \\ J^{t+1} &\infegal& \sum_{i, c_i^t \neq c_i^{t+1}} \; d^2\pa{ X_i, G_{c_i^{t}}^t} + \sum_{i, c_i^t = c_i^{t+1}} \; d^2\pa{ X_i, G_{c_i^{t}}^t} \\ J^{t+1} &\infegal& I^t \end{eqnarray}

Le lemme précédent appliqué à chacune des classes \ensemble{1}{C}, permet d’affirmer que I^{t+1} \infegal J^{t+1}. Par conséquent, la suite \pa{I_t}_{t\supegal 0} est décroissante et minorée par 0, elle est donc convergente.

L’algorithme des centres mobiles cherche à attribuer à chaque point de l’ensemble une classe parmi les C disponibles. La solution trouvée dépend de l’initialisation et n’est pas forcément celle qui minimise l’inertie intra-classe : l’inertie finale est un minimum local. Néanmoins, elle assure que la partition est formée de classes convexes : soit c_1 et c_2 deux classes différentes, on note C_1 et C_2 les enveloppes convexes des points qui constituent ces deux classes, alors \overset{o}{C_1} \cap \overset{o}{C_2} = \emptyset. La figure suivante présente un exemple d’utilisation de l’algorithme des centres mobiles. Des points sont générés aléatoirement dans le plan et répartis en quatre groupes.

../_images/cm.png

C’est une application des centres mobiles avec une classification en quatre classes d’un ensemble aléatoire de points plus dense sur la partie droite du graphe. Les quatre classes ainsi formées sont convexes.

Homogénéité des dimensions#

Les coordonnées des points \left(X_i\right) \in \R^N sont généralement non homogènes : les ordres de grandeurs de chaque dimension sont différents. C’est pourquoi il est conseillé de centrer et normaliser chaque dimension. On note : \forall i \in \intervalle{1}{P}, \; X_i = \vecteur{X_{i,1}}{X_{i,N}} :

\begin{eqnarray*} g_k &=& \pa{EX}_k = \frac{1}{P} \sum_{i=1}^P X_{i,k} \\ v_{kk} &=& \pa{E\left(X-EX\right)^2}_{kk}=\pa{EX^2}_{kk} - g_k^2 \end{eqnarray*}

Les points centrés et normalisés sont :

\forall i \in \intervalle{1}{P}, \; X_i^{\prime}=\left(\dfrac{x_{i,1}-g_{1}}{\sqrt{v_{11}}},...,\dfrac{x_{i,N}-g_{N}}{\sqrt{v_{NN}}}\right)

L’algorithme des centres mobiles est appliqué sur l’ensemble \left( X_{i}^{\prime}\right)_{1\leqslant i\leqslant P}. Il est possible ensuite de décorréler les variables ou d’utiliser une distance dite de Malahanobis définie par d_M\pa{X, Y} = X \, M \, Y' où Y' désigne la transposée de Y et M est une matrice symétrique définie positive. Dans le cas de variables corrélées, la matrice M = \Sigma^{-1} où \Sigma^{-1} est la matrice de variance-covariance des variables aléatoires \pa{X_i}_i.

Améliorations de l’initialisation#

K-means++#

L’article [Arthur2007] montre que l’initialisation aléatoire n’est pas efficace et est sensible aux outliers ou points aberrants. L’étape d’initialisation est remplacée par la suivante :

Algorithme A2 : initialisation k-means++

Cette étape d’initialisation viendra remplacer celle définie dans l’algorithme k-means. On considère un ensemble de points :

X=\left(X_i\right)_{1\leqslant i\leqslant P}\in\left(\R^N\right)^P

A chaque point est associée une classe : \left(c_i\right)_{1\leqslant i\leqslant P}\in\left\{1,...,C\right\}^P.

Pour k centres, on choisit C_1 au hasard dans l’ensemble X. Pour les suivants :

  1. k \leftarrow 2

  2. On choisit aléatoirement G_k \in X avec la probabilité P(x) = \frac{D_{k-1}(x)^2}{\sum_{x\in X}D_{k-1}(x)^2}

  3. k \leftarrow k+1

  4. On revient à l’étape 2 jusqu’à ce que k=C.

La fonction D_k est définie par la distance du point x au centre G_l choisi parmi les k premiers centres. D_k(x) = \min_{1 \infegal l \infegal k} d(x - G_l).

La suite de l’algorithme k-means++ reprend les mêmes étapes que k-means.

Cette initilisation éloigne le prochain centre le plus possibles des centres déjà choisis. L’article montre que :

Théorème T2 : Borne supérieure de l’erreur produite par k-means++

On définit l’inertie par J_(X) = \sum_{i=1}^{P} \; \min_G d^2(X_i, G). Si J_{OPT} définit l’inertie optimale alors \esp{J(X)} \infegal 8 (\ln C + 2) J_{OPT}(X).

La démonstration est disponible dans l’article [Arthur2007].

K-means||#

L’article [Bahmani2012] propose une autre initialisation que K-means++ mais plus rapide et parallélisable.

Algorithme A3 : initialisation k-means||

Cette étape d’initialisation viendra remplacer celle définie dans l’algorithme k-means. On considère un ensemble de points :

X=\left(X_i\right)_{1\leqslant i\leqslant P}\in\left(\R^N\right)^P

A chaque point est associée une classe : \left(c_i\right)_{1\leqslant i\leqslant P}\in\left\{1,...,C\right\}^P.

Pour k centres, on choisit G = \{G_1\} au hasard dans l’ensemble X.

on répète O(\ln D(G, X)) fois :
G' \leftarrow échantillon aléatoire issue de X de probabilité p(x) = l \frac{D(G,x)^2}{\sum_x D(G,x)^2}
G \leftarrow G \cup G'

La fonction D(G,x) est définie par la distance du point x au plus proche centre g \in G : D(g,x) = \min_{g \in G} d(x - g). Cette étape ajoute à l’ensemble des centres G un nombre aléatoire de centres à chaque étape. L’ensemble G contiendra plus de C centres.

  1. Pour tout g \in G, on assigne le poids