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 k1 : 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)
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 k1 : 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 k1 : 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'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}\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 k2 : 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 k2 : 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 k3 : 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 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*}

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

../_images/class_4.png ../_images/class_4_db.png

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 \R^n \longrightarrow \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 : \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 \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.

../_images/herbin1.png ../_images/herbin2.png

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 \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.

../_images/koth1.png ../_images/koth2.png
(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 k4 : 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{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}

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 k5 : 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 k6 : 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.

../_images/class6.png

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.

../_images/liu3.png

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 :

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.

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{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}

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 k7 : 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 \R^2 \longrightarrow \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.

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. 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) PDF, arXiv
[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