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(\mathbb{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(\mathbb{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{split}\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}\end{split}\]

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{\mathbb{R}^N}^P\), \(P\) points de \(\mathbb{R}^N\), le minimum de la quantité \(Q\pa{Y \in \mathbb{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{\mathbb{R}^N}^P\), \(P\) points de \(\mathbb{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\mathbb{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.

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 \mathbb{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(\mathbb{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(\mathbb{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 \(w_g = card \acc{ y | d(x, y) < \min_{h \in G} d(x, h)}\)

2. On clusterise l’ensemble des points \(G\) en \(C\) clusters (avec un k-means classique par exemple)

Au lieu d’ajouter les centres un par un comme dans l’algorithme k-means++, plusieurs sont ajoutés à chaque fois, plus \(l\) est grand, plus ce nombre est grand. Le tirage d’un échantillon aléatoire consiste à inclure chaque point \(x\) avec la probabilité \(p(x) = l \frac{D(G,x)^2}{\sum_x D(G,x)^2}\).

Estimation de probabilités

A partir de cette classification en \(C\) classes, on construit un vecteur de probabilités pour chaque point \(\pa{X_{i}}_{1 \infegal i \infegal P}\) en supposant que la loi de \(X\) sachant sa classe \(c_X\) est une loi normale multidimensionnelle. La classe de \(X_i\) est notée \(c_i\). On peut alors écrire :

\begin{eqnarray*} \forall i \in \intervalle{1}{C}, \; & & \\ G_i &=& E\pa{X \indicatrice{c_X = i}} = \dfrac{\sum_{k=1}^{P} X_k \indicatrice {c_k = i}} {\sum_{k=1}^{P} \indicatrice {c_k = i}} \\ V_i &=& E\pa{XX' \indicatrice{c_X = i}} = \dfrac{\sum_{k=1}^{P} X_k X_k' \indicatrice {c_k = i}} {\sum_{k=1}^{P} \indicatrice {c_k = i}} \\ \pr{c_X = i} &=& \sum_{k=1}^{P} \indicatrice {c_k = i} \\ f\pa{X | c_X = i} &=& \dfrac{1}{\pa{2\pi}^{\frac{N}{2}} \sqrt{\det \pa{V_i}}} \; e^{ - \frac{1}{2} \pa{X - G_i}' \; V_i^{-1} \; \pa{X - G_i} } \\ f\pa{X} &=& \sum_{k=1}^{P} f\pa{X | c_X = i} \pr{c_X = i} \end{eqnarray*}

On en déduit que :

\[\pr{c_X = i |X } = \dfrac{f\pa{X | c_X = i}\pr{c_X = i}} {f\pa{X} }\]

La densité des obervations est alors modélisée par une mélange de lois normales, chacune centrée au barycentre de chaque classe. Ces probabilités peuvent également être apprises par un réseau de neurones classifieur où servir d’initialisation à un algorithme EM.

Sélection du nombre de classes

Critère de qualité

L’algorithme des centres mobiles effectue une classification non supervisée à condition de connaître au préalable le nombre de classes et cette information est rarement disponible. Une alternative consiste à estimer la pertinence des classifications obtenues pour différents nombres de classes, le nombre de classes optimal est celui qui correspond à la classification la plus pertinente. Cette pertinence ne peut être estimée de manière unique, elle dépend des hypothèses faites sur les éléments à classer, notamment sur la forme des classes qui peuvent être convexes ou pas, être modélisées par des lois normales multidimensionnelles, à matrice de covariances diagonales, … Les deux critères qui suivent sont adaptés à l’algorithme des centres mobiles. Le critère de Davies-Bouldin (voir [Davies1979]) est minimum lorsque le nombre de classes est optimal.

\begin{eqnarray} DB &=& \dfrac{1}{C} \; \sum_{i=1}^{C} \; \max_{i \neq j} \; \dfrac{\sigma_i + \sigma_j}{ d\pa{C_i,C_j}} \end{eqnarray}

Avec :

\(C\)	nombre de classes
\(\sigma_i\)	écart-type des distances des observations de la classe \(i\)
\(C_i\)	centre de la classe \(i\)

Le critère de Goodman-Kruskal (voir [Goodman1954]) est quant à lui maximum lorsque le nombre de classes est optimal. Il est toutefois plus coûteux à calculer.

\begin{eqnarray} GK &=& \dfrac{S^+ - S^-} { S^+ + S^-} \end{eqnarray}

Avec :

\begin{eqnarray*} S^+ &=& \acc{ \pa{q,r,s,t} \sac d\pa{q,r} < d\pa{s,t} } \\ S^- &=& \acc{ \pa{q,r,s,t} \sac d\pa{q,r} < d\pa{s,t} } \end{eqnarray*}

Où \(\pa{q,r}\) sont dans la même classe et \(\pa{s,t}\) sont dans des classes différentes.

Classification en quatre classes : nombre de classes sélectionnées par le critère de Davies-Bouldin dont les valeurs sont illustrées par le graphe apposé à droite.

Maxima de la fonction densité

L’article [Herbin2001] propose une méthode différente pour estimer le nombre de classes, il s’agit tout d’abord d’estimer la fonction densité du nuage de points qui est une fonction de \(\mathbb{R}^n \longrightarrow \mathbb{R}\). Cette estimation est effectuée au moyen d’une méthode non paramètrique telle que les estimateurs à noyau (voir [Silverman1986]) Soit \(\vecteur{X_1}{X_N}\) un nuage de points inclus dans une image, on cherche à estimer la densité \(f_H\pa{x}\) au pixel \(x\) :

\[\hat{f}_H\pa{x} = \dfrac{1}{N} \; \sum_{i=1}^{N} \; \dfrac{1}{\det H} \; K\pa{ H^{-1} \pa{x - X_i}}\]

Où :

\[K\pa{x} = \dfrac{1}{ \pa{2 \pi}^{ \frac{d}{2}} } \; e^{ - \frac{ \norme{x}^2 } {2} }\]

\(H\) est un paramètre estimée avec la règle de Silverman. L’exemple utilisé dans cet article est un problème de segmentation d’image qui ne peut pas être résolu par la méthode des nuées dynamiques puisque la forme des classes n’est pas convexe, ainsi que le montre la figure suivante. La fonction de densité \(f\) est seuillée de manière à obtenir une fonction \(g : \mathbb{R}^n \longrightarrow \acc{0,1}\) définie par :

\[g \pa{x} = \indicatrice{f\pa{x} \supegal s}\]

L’ensemble \(g^{-1}\pa{\acc{1}} \subset \mathbb{R}^n\) est composée de \(N\) composantes connexes notées \(\vecteur{C_1}{C_N}\), la classe d’un point \(x\) est alors l’indice de la composante connexe à la laquelle il appartient ou la plus proche le cas échéant.

Exemple de classification non supervisée appliquée à un problème de segmentation d’image, la première figure montre la densité obtenue, la seconde figure illustre la classification obtenue, figure extraite de [Herbin2001]. Cette méthode paraît néanmoins difficilement applicable lorsque la dimension de l’espace vectoriel atteint de grande valeur. L’exemple de l’image est pratique, elle est déjà découpée en région représentées par les pixels, l’ensemble \(g^{-1}\pa{\acc{1}}\) correspond à l’ensemble des pixels \(x\) pour lesquels \(f\pa{x} \supegal s\).

Décroissance du nombre de classes

L’article [Kothari1999] propose une méthode permettant de faire décroître le nombre de classes afin de choisir le nombre approprié. L’algorithme des centres mobiles proposent de faire décroître l’inertie notée \(I\) définie pour un ensemble de points noté \(X = \vecteur{x_1}{x_N}\) et \(K\) classes. La classe d’un élément \(x\) est notée \(C\pa{x}\). Les centres des classes sont notés \(Y = \vecteur{y_1}{y_K}\). L’inertie de ce nuage de points est définie par :

\[I = \sum_{x \in X} \; \norme{ x - y_{C\pa{x} }}^2\]

On définit tout d’abord une distance \(\alpha \in \mathbb{R}^+\), puis l’ensemble \(V\pa{y,\alpha} = \acc{ z \in Y \sac d\pa{y,z} \infegal \alpha }\), \(V\pa{y,\alpha}\) est donc l’ensemble des voisins des centres dont la distance avec \(y\) est inférieur à \(\alpha\). L’article [Kothari1999] propose de minimiser le coût \(J\pa{\alpha}\) suivant :

\[J\pa{\alpha} = \sum_{x \in X} \; \norme{ x - y_{C\pa{x} }}^2 + \sum_{x \in X} \; \sum_{y \in V\pa{y_{C\pa{x}}, \alpha} } \; \lambda\pa{y} \, \norme{ y - y_{C\pa{x}}}^2\]

Lorsque \(\alpha\) est nul, ce facteur est égal à l’inertie : \(I = J\pa{0}\) et ce terme est minimal lorsqu’il y a autant de classes que d’éléments dans \(X\). Lorsque \(\alpha\) tend vers l’infini, \(J\pa{\alpha} \rightarrow J\pa{\infty}\) où :

\[J\pa{\infty} = \sum_{x \in X} \; \norme{ x - y_{C\pa{x} }}^2 + \sum_{x \in X} \; \sum_{y \in Y} \; \lambda\pa{y} \, \norme{ y - y_{C\pa{x}}} ^2\]

Ici encore, il est possible de montrer que ce terme \(J\pa{\infty}\) est minimal lorsqu’il n’existe plus qu’une seule classe. Le principe de cette méthode consiste à faire varier le paramètre \(\alpha\), plus le paramètre \(\alpha\) augmente, plus le nombre de classes devra être réduit. Néanmoins, il existe des intervalles pour lequel ce nombre de classes est stable, le véritable nombre de classes de l’ensemble \(X\) sera considéré comme celui correspondant au plus grand intervalle stable.

(a)	(b)

Evolutation du nombre de classes en fonction du paramètre \(\alpha\) lors de la minimisation du critère \(J\pa{\alpha}\), figure extraite de [Kothari1999]. La première image représente le nuage de points illustrant quatre classes sans recouvrement. La seconde image montre que quatre classes est l’état le plus longtemps stable lorsque \(\alpha\) croît.

Le coût \(J\pa{\alpha}\) est une somme de coût dont l’importance de l’un par rapport à l’autre est contrôle par les paramètres \(\lambda\pa{y}\). Le problème de minimisation de \(J\pa{\alpha}\) est résolu par l’algorithme qui suit. Il s’appuie sur la méthode des multiplicateurs de Lagrange.

Algorithme A4 : sélection du nombre de classes

(voir [Kothari1999]) Les notations sont celles utilisés dans les paragraphes précédents. On suppose que le paramètre \(\alpha\) évolue dans l’intervalle \(\cro{\alpha_1, \alpha_2}\) à intervalle régulier \(\alpha_t\). Le nombre initial de classes est noté \(K\) et il est supposé surestimer le véritable nombre de classes. Soit \(\eta \in \left]0,1\right[\), ce paramètre doit être choisi de telle sorte que dans l’algorithme qui suit, l’évolution des centres \(y_k\) soit autant assurée par le premier de la fonction de coût que par le second.

initialisation

\[\alpha \longleftarrow \alpha_1\]

On tire aléatoirement les centres des \(K\) classes \(\vecteur{y_1}{y_K}\).

préparation

On définit les deux suites entières \(\vecteur{c^1_1}{c^1_K}\), \(\vecteur{c^2_1}{c^2_K}\), et les deux suites de vecteur \(\vecteur{z^1_1}{z^1_K}\), \(\vecteur{z^2_1}{z^2_K}\).

\[\begin{split}\begin{array}{rlll} \forall k, & c^1_k &=& 0 \\ \forall k, & c^2_k &=& 0 \\ \forall k, & z^1_k &=& 0 \\ \forall k, & z^2_k &=& 0 \end{array}\end{split}\]

calcul des mises à jour

for i in \(1..N\)

Mise à jour d’après le premier terme de la fonction de coût \(J\pa{\alpha}\).

\(w \longleftarrow \underset{1 \infegal l \infegal K}{\arg \min} \; \norme{x_i - y_l}^2\)

\(z^1_w \longleftarrow z^1_w + \eta \pa{ x_i - y_w}\)

\(c^1_w \longleftarrow c^1_w + 1\)

Mise à jour d’après le second terme de la fonction de coût \(J\pa{\alpha}\)

for v in \(1..k\)

if \(\norme{y_v - y_w} < \alpha\)

\(z^2_v \longleftarrow z^2_v - \pa{ y_v - y_w}\)

\(c^2_v \longleftarrow c^2_v + 1\)

for v in \(1..k\)

\(\lambda_v \longleftarrow \frac{ c^2_v \norme{z^1_v} } { c^1_v \norme{z^2_v} }\)

\(y_v \longleftarrow y_v + z^1_v + \lambda_v z^2_v\)

convergence

Tant que l’étape précédente n’a pas convergé vers une version stable des centres, \(y_k\), retour à l’étape précédente. Sinon, tous les couples de classes \(\pa{i,j}\) vérifiant \(\norme{y_i - y_j} > \alpha\) sont fusionnés : \(\alpha \longleftarrow \alpha + \alpha_t\). Si \(\alpha \infegal \alpha2\), retour à l’étape de préparation.

terminaison

Le nombre de classes est celui ayant prévalu pour le plus grand nombre de valeur de \(\alpha\).

Extension des nuées dynamiques

Classes elliptiques

La version de l’algorithme des nuées dynamique proposée dans l’article [Cheung2003] suppose que les classes ne sont plus de forme circulaire mais suivent une loi normale quelconque. La loi de l’échantillon constituant le nuage de points est de la forme :

\[f\pa{x} = \sum_{i=1}^{N} \; p_i \; \dfrac{1}{\pa{2 \pi}^{\frac{d}{2}}\sqrt{\det \Sigma_i}} \; exp \pa{-\frac{1}{2} \pa{x-\mu_i}' \Sigma_i^{-1} \pa{x-\mu_i} }\]

Avec \(sum_{i=1}^{N} \; p_i = 1\). On définit :

\[G\pa{x, \mu, \Sigma} = \dfrac{1}{\pa{2 \pi}^{\frac{d}{2}}\sqrt{\det \Sigma}} \; exp \pa{-\frac{1}{2} \pa{x-\mu}' \Sigma^{-1} \pa{x-\mu} }\]

L’algorithme qui suit a pour objectif de minimiser la quantité pour un échantillon \(\vecteur{X_1}{X_K}\) :

\[I = \sum_{i=1}^{N}\sum_{k=1}^{K} \indicatrice{ i = \underset{1 \infegal j \infegal N}{\arg \max} G\pa{X_k, \mu_j,\Sigma_j} } \; \ln \cro{ p_i G\pa{ X_k, \mu_i, \Sigma_i } }\]

Algorithme A5 : nuées dynamiques généralisées

Les notations sont celles utilisées dans ce paragraphe. Soient \(\eta\), \(\eta_s\) deux réels tels que \(\eta > \eta_s\). La règle préconisée par l’article [Cheung2003] est \(\eta_s \sim \frac{\eta}{10}\).

initialisation

\(t \longleftarrow 0\). Les paramètres \(\acc{p_i^0, \mu_i^0, \Sigma_i^0 \sac 1 \infegal i \infegal N}\) sont initialisés grâce à un algorithme des k-means ou FSCL. \(\forall i, \; p_i^0 = \frac{1}{N}\) et \(\beta_i^0 = 0\).

récurrence

Soit \(X_k\) choisi aléatoirement dans \(\vecteur{X_1}{X_K}\).

\[i = \underset{1 \infegal i \infegal N}{\arg \min} \; G\pa{X_k, \mu_i^t, \Sigma_i^t}\]

for i in \(1..N\)

\(\mu_i^{t+1} = \mu_i^t + \eta \, \pa{\Sigma_i^t}^{-1} \, \pa{ X_k - \mu_i^t}\)

\(\beta_i^{t+1} = \beta_i^t + \eta \, \pa{1 - \alpha_i^t}\)

\(\Sigma^{t+1}_i = \pa{1 - \eta_s} \, \Sigma_i^t + \eta_s \, \pa{ X_k - \mu_i^t} \pa{ X_k - \mu_i^t}'\)

for i in \(1..N\)

\(p^{t+1}_i = \frac{ e^{ \beta_i^{t+1} } } { \sum_{j=1}^{N} e^{ \beta_j^{t+1} } }\)

\(t \longleftarrow t + 1\)

terminaison

Tant que \(\underset{1 \infegal i \infegal N}{\arg \min} \; G\pa{X_k, \mu_i^t, \Sigma_i^t}\) change pour au moins un des points \(X_k\).

Lors de la mise à jour de \(\Sigma^{-1}\), l’algorithme précédent propose la mise à jour de \(\Sigma_i\) alors que le calcul de \(G\pa{., \mu_i, \Sigma_i}\) implique \(\Sigma_i^{-1}\), par conséquent, il est préférable de mettre à jour directement la matrice \(\Sigma^{-1}\) :

\[\pa{\Sigma^{t+1}_i}^{-1} = \frac{ \pa{\Sigma_i^t}^{-1} } {1 - \eta_s} \cro{I - \frac{ \eta_s \pa{ X_k - \mu_i^t} \pa{ X_k - \mu_i^t}' \pa{\Sigma_i^t}^{-1} } {1 - \eta_s + \eta_s \pa{ X_k - \mu_i^t}' \, \pa{\Sigma_i^t}^{-1}\pa{ X_k - \mu_i^t} } }\]

Rival Penalized Competitive Learning (RPCL)

L’algorithme suivant développé dans [Xu1993], est une variante de celui des centres mobiles. Il entreprend à la fois la classification et la sélection du nombre optimal de classes à condition qu’il soit inférieur à une valeur maximale à déterminer au départ de l’algorithme. Un mécanisme permet d’éloigner les centres des classes peu pertinentes de sorte qu’aucun point ne leur sera affecté.

Algorithme A6 : RPCL

Soient \(\vecteur{X_1}{X_N}\), \(N\) vecteurs à classer en au plus \(T\) classes de centres \(\vecteur{C_1}{C_T}\). Soient deux réels \(\alpha_r\) et \(\alpha_c\) tels que \(0 < \alpha_r \ll \alpha_c < 1\).

initialisation

Tirer aléatoirement les centres \(\vecteur{C_1}{C_T}\).

for j in \(1..C\)

\(n_j^0 \longleftarrow 1\)

calcul de poids

Choisir aléatoirement un point \(X_i\).

for j in \(1..C\)

\(\gamma_j = \dfrac{n_j}{ \sum_{k=1}^{C} n_k}\)

for j in \(1..C\)

\(u_j =\)

1 si \(j \in \underset{k}{\arg \min} \; \cro {\gamma_k \; d\pa{X_i,C_k} }\)

-1 si \(j \in \underset{j \neq k}{\arg \min} \; \cro {\gamma_k \; d\pa{X_i,C_k} }\)

0 sinon

mise à jour

for j in \(1..C\)

\(C_j^{t+1} \longleftarrow C_j^t + \left \{ \begin{array}{ll} \alpha_c \pa{X_i - C_j} & \text{si } u_j = 1 \\ - \alpha_r \pa{X_i - C_j} & \text{si } u_j = -1 \\ 0 & \text{sinon} \end{array} \right.\)

\(n_j^t + \left \{ \begin{array}{ll} 1 & \text{si } u_j = 1 \\ 0 & \text{sinon} \end{array} \right.\)

\(t \longleftarrow t+1\)

terminaison

S’il existe un indice \(j\) pour lequel \(C^{t+1}_j \neq C^t_j\) alors retourner à l’étape de calcul de poids ou que les centres des classes jugées inutiles ont été repoussés vers l’infini.

Pour chaque point, le centre de la classe la plus proche en est rapproché tandis que le centre de la seconde classe la plus proche en est éloigné mais d’une façon moins importante (condition \(\alpha_r \ll \alpha_c\)). Après convergence, les centres des classes inutiles ou non pertinentes seront repoussés vers l’infini. Par conséquent, aucun point n’y sera rattaché.

L’algorithme doit être lancé plusieurs fois. L’algorithme RPCL peut terminer sur un résultat comme celui de la figure suivante où un centre reste coincé entre plusieurs autres. Ce problème est moins important lorsque la dimension de l’espace est plus grande.

Application de l’algorithme RPCL : la classe 0 est incrusté entre les quatre autres et son centre ne peut se « faufiler » vers l’infini.

RPCL-based local PCA

L’article [Liu2003] propose une extension de l’algorithme RPCL et suppose que les classes ne sont plus de forme circulaire mais suivent une loi normale quelconque. Cette méthode est utilisée pour la détection de ligne considérées ici comme des lois normales dégénérées en deux dimensions, la matrice de covariance définit une ellipse dont le grand axe est très supérieur au petit axe, ce que montre la figure suivante. Cette méthode est aussi présentée comme un possible algorithme de squelettisation.

Figure extraite de [Liu2003], l’algorithme est utilisé pour la détection de lignes considérées ici comme des lois normales dont la matrice de covariance définit une ellipse dégénérée dont le petit axe est très inférieur au grand axe. Les traits fin grisés correspondent aux classes isolées par l’algorithme RPCL-based local PCA.

On modélise le nuage de points par une mélange de lois normales :

\[f\pa{x} = \sum_{i=1}^{N} \; p_i \; \dfrac{1}{\pa{2 \pi}^{\frac{d}{2}}\sqrt{\det \Sigma_i}} \; exp \pa{-\frac{1}{2} \pa{x-\mu_i}' \Sigma_i^{-1} \pa{x-\mu_i} }\]

Avec \(\sum_{i=1}^{N} \; p_i = 1\).

On suppose que le nombre de classes initiales \(N\) surestime le véritable nombre de classes. L’article [Liu2003] s’intéresse au cas particulier où les matrices de covariances vérifient \(\Sigma_i = \zeta_i \, I + \sigma_i \, \phi_i \phi_i'\) avec \(\zeta_i > 0, \; \sigma_i > 0, \; \phi_i' \phi_i = 1\).

On définit également :

\[G\pa{x, \mu, \Sigma} = \dfrac{1}{\pa{2 \pi}^{\frac{d}{2}}\sqrt{\det \Sigma}} \; exp \pa{-\frac{1}{2} \pa{x-\mu}' \Sigma^{-1} \pa{x-\mu} }\]

L’algorithme utilisé est similaire à l’algortihme RPCL. La distance \(d\) utilisée lors de l’étape de calcul des poids afin de trouver la classe la plus probable pour un point donné \(X_k\) est remplacée par l’expression :

\[d\pa{X_k, classe \, i} = - \ln { p_i^t \, G\pa{X_k, \, \mu_i^t, \, \Sigma^t_i } }\]

L’étape de mise à jour des coefficients est remplacée par :

\[\begin{split}x^{t+1} \longleftarrow x^t + \left \{ \begin{array}{ll} \alpha_c \nabla x^t & \text{si } u_j = 1 \\ - \alpha_r \nabla x^t & \text{si } u_j = -1 \\ 0 & \text{sinon} \end{array} \right.\end{split}\]

Où \(x^t\) joue le rôle d’un paramètre et est remplacé successivement par \(p_i^t\), \(\mu_i^t\), \(\zeta_i^t\), \(\sigma^t_i\), \(\phi^t_i\) :

\[\begin{split}\begin{array}{lll} \nabla p_i^t &=& - \frac{1}{p_i^t} \\ \nabla \mu_i^t &=& - \pa{ X_k - \mu_i^t} \\ \nabla \zeta_i^t &=& \frac{1}{2} \; tr\cro{ \pa{\Sigma_i^t}^{-1} \, \pa{ I - \pa{ X_k - \mu_i^t} \pa{ X_k - \mu_i^t}' \pa{\Sigma_i^t}^{-1} } } \\ \nabla \sigma_i^t &=& \frac{1}{2} \; \pa{\phi_i^t}' \pa{\Sigma_i^t}^{-1} \pa{ I - \pa{ X_k - \mu_i^t} \pa{ X_k - \mu_i^t}' \pa{\Sigma_i^t}^{-1} } \phi_i^t \\ \nabla \phi_i^t &=& \sigma_i^t \pa{\Sigma_i^t}^{-1} \pa{ I - \pa{ X_k - \mu_i^t} \pa{ X_k - \mu_i^t}' \pa{\Sigma_i^t}^{-1} } \phi_i^t \\ \end{array}\end{split}\]

Frequency Sensitive Competitive Learning (FSCL)

L’algorithme Frequency Sensitive Competitive Learning est présenté dans [Balakrishnan1996]. Par rapport à l’algorithme des centres mobiles classique, lors de l’estimation des centres des classes, l’algorithme évite la formation de classes sous-représentées.

Algorithme A7 : FSCL

Soit un nuage de points \(\vecteur{X_1}{X_N}\), soit \(C\) vecteurs \(\vecteur{\omega_1}{\omega_C}\) initialisés de manière aléatoires. Soit \(F : \pa{u,t} \in \mathbb{R}^2 \longrightarrow \mathbb{R}^+\) croissante par rapport à \(u\). Soit une suite de réels \(\vecteur{u_1}{u_C}\), soit une suite \(\epsilon\pa{t} \in \cro{0,1}\) décroissante où \(t\) représente le nombre d’itérations. Au début \(t \leftarrow 0\).

meilleur candidat

Pour un vecteur \(X_k\) choisi aléatoirement dans l’ensemble \(\vecteur{X_1}{X_N}\), on détermine :

\[i^* \in \arg \min \acc{ D_i = F\pa{u_i,t} \, d\pa{X_k, \omega_i} }\]

mise à jour

\(\omega_{i^*} \pa{t+1} \longleftarrow \omega_{i^*} \pa{t} + \epsilon\pa{t} \pa { X_k - \omega_{i^*} \pa{t} }\)

\(t \longleftarrow t+1\)

\(u_{i^*} \longleftarrow u_{i^*} + 1\)

Retour à l’étape précédente jusqu’à ce que les nombres \(\frac{u_i}{\sum_{i}u_i}\) convergent.

Exemple de fonctions pour \(F\), \(\epsilon\) (voir [Balakrishnan1996]) :

\begin{eqnarray*} F\pa{u,t} &=& u \, \beta e^{-t/T} \text{ avec } \beta = 0,06 \text{ et } 1/T = 0,00005 \\ \epsilon\pa{t} &=& \beta \, e^{ - \gamma t } \text{ avec } \gamma = 0,05 \end{eqnarray*}

Cet algorithme ressemble à celui des cartes topographiques de Kohonen sans toutefois utiliser un maillage entre les neurones (ici les vecteurs \(\omega_i\)). Contrairement à l’algorithme RPCL, les neurones ne sont pas repoussés s’ils ne sont pas choisis mais la fonction croissante \(F\pa{u,t}\) par rapport à \(u\) assure que plus un neurone est sélectionné, moins il a de chance de l’être, bien que cet avantage disparaisse au fur et à mesure des itérations.

k-means norme L1

L’algorithme dans sa version la plus courante optimise l’inertie définie par \(\sum_{i=1}^P \; d^2\left(X_i, G_{c_i^t}^t\right)\), qui est en quelque sorte une inertie L2. Que devriendrait l’algorithme si la norme choisie était une norme L1, il faudrait alors choisir à chaque itération t des points qui minimise la quantité : \(\sum_{i=1}^P \; d_1\left(X_i, G_{c_i^t}^t\right)\) où \(d_1\) est la norme L1 entre deux points X,Y : \(d_1(X, Y) = \sum_i |X_i - Y_i|\). Avant de continuer, on rappelle un théorème :

propriété P1 : Médiane et valeur absolue

Soit \(A=(x_1, ..., x_n)\) un ensembl de n réels quelconque. On note \(m=med(x_1, ..., x_n)\) la médiane de l’ensemble de points A. Alors la médiane m minimise la quantité \(\sum_{i=1}^n |m-x_i|\).

C’est cette propriété qui est utilisée pour définir ce qu’est la régression quantile et sa démonstration est présentée à la page Médiane et valeur absolue. Il ne reste plus qu’à se servir de ce résultat pour mettre à jour l’algorithme centre mobile, k-means. L’étape qui consiste à affecter un point à un cluster représenté par un point ne pose pas de problème si on utilise cette nouvelle norme. Il ne reste plus qu’à déterminer le point qui représente un cluster sachant les points qui le constituent. Autrement dit, il faut déterminer le point qui minimiser la pseudo-inertie définie comme suit pour un ensemble de points \((X_1, ..., X_n)\) appartenant à un espace vectoriel de dimension k.

\[I(G,X_1,...,X_n) = \norm{G - X_i}_1 = \sum_{i=1}^n \sum_{k=1}^d \abs{G_k - X_{ik}}\]

On cherche le point G qui minimise la quantité \(I(G,X_1,...,X_n)\). Comme \(\sum_{i=1}^n \sum_{k=1}^d \abs{G_k - X_{ik}} = \sum_{k=1}^d \sum_{i=1}^n \abs{G_k - X_{ik}}\), on en déduit qu’on peut chercher la coordonnée \(G_k\) indépendemment les unes des autres. On en déduit que le barycentre de norme L1 d’un ensemble de points dans un espace vectoriel de dimension d a pour coordonnées les d médianes extraites sur chacune des dimensions. L’algorithme est implémenté dans le module mlinsights en s’inspirant du code KMeans.

Bibliographie

[Arthur2007] (1,2)

k-means++: the advantages of careful seeding (2007), Arthur, D.; Vassilvitskii, S., Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms. Society for Industrial and Applied Mathematics Philadelphia, PA, USA. pp. 1027–1035. 2006-13.pdf.

[Balakrishnan1996] (1,2)

Comparative performance of the FSCL neural net and K-means algorithm for market segmentation (1996), P. V. Sundar Balakrishnan, Martha Cooper, Varghese S. Jacob, Phillip A. Lewis, European Journal of Operation Research, volume 93, pages 346-357

[Bahmani2012]

Scalable K-Means++ (2012), Bahman Bahmani, Benjamin Moseley, Andrea Vattani, Ravi Kumar, Sergei Vassilvitskii, Proceedings of the VLDB Endowment (PVLDB), Vol. 5, No. 7, pp. 622-633 (2012) vldb12-kmpar.pdf, arXiv.1203.6402

[Cheung2003] (1,2)

\(k^*\)-Means: A new generalized k-means clustering algorithm (2003), Yiu-Ming Cheung, Pattern Recognition Letters, volume 24, 2883-2893

[Davies1979]

A cluster Separation Measure (1979), D. L. Davies, D. W. Bouldin, IEEE Trans. Pattern Analysis and Machine Intelligence (PAMI), volume 1(2)

[Goodman1954]

Measures of associations for cross-validations (1954), L. Goodman, W. Kruskal, J. Am. Stat. Assoc., volume 49, pages 732-764

[Herbin2001] (1,2)

Estimation of the number of clusters and influence zones (2001), M. Herbin, N. Bonnet, P. Vautrot, Pattern Recognition Letters, volume 22, pages 1557-1568

[Kothari1999] (1,2,3,4)

On finding the number of clusters (1999), Ravi Kothari, Dax Pitts, Pattern Recognition Letters, volume 20, pages 405-416

[Liu2003] (1,2,3)

Strip line detection and thinning by RPCL-based local PCA (2003), Zhi-Yong Liu, Kai-Chun Chiu, Lei Xu, Pattern Recognition Letters volume 24, pages 2335-2344

[Silverman1986]

Density Estimation for Statistics and Data Analysis (1986), B. W. Silverman, Monographs on Statistics and Applied Probability, Chapman and Hall, London, volume 26

[Xu1993]

Rival penalized competitive learning for clustering analysis, rbf net and curve detection (1993), L. Xu, A. Krzyzak, E. Oja, IEEE Trans. Neural Networks, volume (4), pages 636-649