Cours - Classification automatique

Ce chapitre correspond à 2 séances de cours : la première porte sur K-means et sur la Classification ascendante hiérarchique, la seconde sur la Classification spectrale.

[Diapositives du cours]

Généralités

La classification automatique (cluster analysis ou clustering en anglais) cherche à répartir un ensemble donné de \(N\) observations en groupes (catégories, classes, taxons, clusters) de façon à regrouper les observations similaires et à séparer les observations dissimilaires. Dans la mesure où les notions de similarité et de groupe peuvent être explicitées de multiples façons, de nombreuses méthodes de classification automatique ont été proposées depuis les années 1930. Des synthèses partielles peuvent être trouvées, par exemple, dans [JMF99] ou [BBS14]. Le choix d’une méthode adéquate exige une connaissance des données et de la nature des groupes recherchés.

En général, nous pouvons parler de classification automatique si aucune information n’est disponible concernant l’appartenance de certaines données à certaines classes connues a priori. Par ailleurs, le nombre de groupes recherchés peut être connu a priori ou non.

Une première distinction entre méthodes de classification automatique peut être faite suivant leur objectif. La plupart des méthodes visent à obtenir un partitionnement des données, comme dans l’exemple suivant, plus éventuellement à trouver une donnée « représentative » (« prototype ») par groupe :

Exemple de partitionnement de données bidimensionnelles en 3 groupes

Fig. 60 Exemple de partitionnement de données bidimensionnelles en 3 groupes

D’autres méthodes cherchent plutôt à obtenir une hiérarchie de regroupements, qui fournit une information plus riche concernant la structure de similarité des données. Noter qu’à partir d’une telle hiérarchie il est facile d’extraire plusieurs partitionnements, à des niveaux de « granularité » différents, comme dans l’exemple suivant :

Exemple de regroupement hiérarchique

Fig. 61 Exemple de regroupement hiérarchique qui peut servir à obtenir plusieurs partitionnements

Dans les deux cas, la classification automatique permet de

  • mettre en évidence une structure simple (partitionnement ou hiérarchie de groupes) dans un ensemble de données ;

  • résumer un ensemble de données par les représentants des groupes (« prototypes »), éventuellement à plusieurs niveaux hiérarchiques.

Une seconde distinction entre méthodes de classification automatique peut être faite suivant la nature des données à regrouper : données numériques, données catégorielles ou données mixtes. Si la plupart des méthodes s’intéressent aux données numériques, certaines peuvent traiter directement des données catégorielles. Toutefois, les données catégorielles sont fréquemment transformées en données numériques (par ex. par une analyse factorielle des correspondances multiples) avant d’être traitées par des méthodes de classification adaptées aux données numériques.

Une autre distinction correspond à la représentation des données. De nombreuses méthodes de classification automatique travaillent sur des représentations vectorielles qui, en plus d’accepter différentes métriques, permettent de définir des centres de gravité de groupes, des densités de probabilité, des intervalles, des sous-espaces, etc. Ces méthodes sont en général basées sur l’utilisation de prototypes et leur complexité est en général de \(O(N)\). Certaines méthodes se satisfont, en revanche, d’une simple structure métrique (distances) sur l’espace des données ; cela leur confère une grande généralité mais également une complexité plus élevée, \(\geq O(N^2)\).

Exemple de groupes non exclusifs

Fig. 62 A quel groupe appartiennent les données entourées ?

Une distinction importante concerne la nature des groupes recherchés : sont-il mutuellement exclusifs ou non ? Sont-ils « nets » (crisp, une observation appartient ou n’appartient pas à un groupe) ou « flous » (fuzzy, une observation peut appartenir à différents degrés à plusieurs groupes) ? Si la plupart des méthodes de classification automatique s’intéressent aux partitionnement exclusifs, certaines permettent d’obtenir des groupes non exclusifs. Aussi, des extensions floues existent pour des méthodes de classification automatiques bien connues (par ex. fuzzy c-means pour K-means, voir [JMF99]). Les extensions floues peuvent avoir une convergence plus robuste vers une solution que les méthodes « nettes » de départ. Il est possible de passer des groupes flous à des groupes nets en affectant chaque donnée (ou observation) au groupe auquel elle appartient le plus.

Enfin, une autre distinction importante peut être faite suivant le critère de regroupement. Certaines méthodes cherchent des groupes « compacts » et relativement éloignés entre eux, comme dans cet exemple de partitionnement. D’autres méthodes s’intéressent à des groupes denses (et non nécessairement « compacts ») séparés par des régions moins denses, comme dans l’exemple suivant :

Groupes denses séparés par des régions moins denses

Fig. 63 Exemple de groupes denses séparés par des régions moins denses

Il faut noter que ce critère de regroupement n’est pas toujours explicite, pourtant son impact sur les résultats obtenus est majeur.

K-means

Considérons un ensemble \(\mathcal{E}\) de \(N\) données décrites par \(p\) variables à valeurs dans \(\mathbb{R}\) et \(d\) une distance sur \(\mathbb{R}^p\). Pour regrouper ces données en \(k\) groupes disjoints \(\mathcal{E}_1,\ldots,\mathcal{E}_k\), inconnus a priori, on s’intéresse souvent à un critère qui correspond à la somme des inerties intra-classe des groupes :

(21)\[\phi_\mathcal{E}(\mathcal{C}) = \sum_{j=1}^k \sum_{\mathbf{x}_i \in \mathcal{E}_j} d^2(\mathbf{x}_i,\mathbf{m}_j)\]

avec \(\mathcal{C} = \{\mathbf{m}_j, 1 \leq j \leq k\}\) l’ensemble des centres de gravité des \(k\) groupes. Pour \(\mathcal{E}\) et \(k\) donnés, plus la valeur de \(\phi_\mathcal{E}(\mathcal{C})\) est faible, plus les groupes sont « compacts » autour de leurs centres et donc meilleure est la qualité du partitionnement obtenu.

La somme des inerties intra-classe peut s’écrire aussi

\[\phi_\mathcal{E}(\mathcal{C}) = \sum_{1 \leq l \leq N} d^2(\mathbf{x}_l,\mathbf{m}_{C(l)})\]

\({C(l)}\) est l’indice du groupe dont fait partie \(\mathbf{x}_l\). Cela permet de voir que, pour un ensemble \(\mathcal{E}\) donné, les variables dont dépend la valeur du critère sont les positions des centres, \(\{\mathbf{m}_j, 1 \leq j \leq k\}\).

Trouver le minimum global de la fonction (21) est un problème NP-difficile, mais on dispose d’algorithmes de complexité polynomiale dans le nombre de données \(N\) qui produisent une solution en général sous-optimale.

Un tel algorithme est l’algorithme des centres mobiles décrit ci-dessous :

Entrées : ensemble \(\mathcal{E}\) de \(N\) données de \(\mathbb{R}^p\) ;

Sorties : \(k\) groupes (clusters) disjoints \(\mathcal{E}_1,\ldots,\mathcal{E}_k\) et ensemble \(\mathcal{C}\) de leurs centres ;

  1. Initialisation aléatoire des centres \(\mathbf{m}_j\), \(1 \leq j \leq k\) ;

  2. tant que (centres non stabilisés) faire

    1. Affectation de chaque donnée au groupe du centre le plus proche ;

    2. Remplacement des centres par les centres de gravité des groupes ;

fin tant que

Des étapes successives dans l’application de cet algorithme à des données issues de 7 lois normales bidimensionnelles peuvent être visualisées dans la figure suivante. La distance euclidienne classique est employée et la classification est réalisée avec avec \(k=7\) centres, dont la position initiale (choisie aléatoirement) est indiquée dans la première image.

Étapes successives dans l'application de *K-means*

Fig. 64 Étapes successives dans l’application de K-means

Dans cet exemple, la solution obtenue est optimale et est atteinte après seulement 2 itérations.

Question : La valeur minimale que peut atteindre \(\phi_\mathcal{E}(\mathcal{C})\)

  1. diminue avec l’augmentation de \(k\),

  2. ne dépend pas de \(k\),

  3. augmente avec l’augmentation de \(k\).


On peut montrer que la valeur de \(\phi_\mathcal{E}(\mathcal{C})\) diminue lors de chacune des deux étapes du processus itératif (affectation de chaque donnée à un groupe, recalcul des centres) :

  1. Affectation de chaque donnée au groupe du centre le plus proche : \(\mathbf{x}_i\) passe du groupe de centre \(\mathbf{m}_p\) au groupe de centre \(\mathbf{m}_q\) si \(d^2(\mathbf{x}_i,\mathbf{m}_p) > d^2(\mathbf{x}_i,\mathbf{m}_q)\), donc

    \[d^2(\mathbf{x}_i,\mathbf{m}_p) + \sum_{l \neq i} d^2(\mathbf{x}_l,\mathbf{m}_{C(l)}) > d^2(\mathbf{x}_i,\mathbf{m}_q) + \sum_{l \neq i} d^2(\mathbf{x}_l,\mathbf{m}_{C(l)})\]
  2. Remplacement des anciens centres par les centres de gravité des groupes : si \(\widetilde{\mathbf{m}}_j\) est l’ancien centre du groupe \(j\) et \(\mathbf{m}_j\) le nouveau, alors

    \[\begin{split}\begin{aligned} \sum_{\mathbf{x}_i \in \mathcal{E}_j} d^2(\mathbf{x}_i,\widetilde{\mathbf{m}}_j) & = \sum_{\mathbf{x}_i \in \mathcal{E}_j} \left\|\mathbf{x}_i - \mathbf{m}_j + \mathbf{m}_j - \widetilde{\mathbf{m}}_j \right\|^2 \\ & = \sum_{\mathbf{x}_i \in \mathcal{E}_j} \left\|\mathbf{x}_i - \mathbf{m}_j \right\|^2 + \sum_{\mathbf{x}_i \in \mathcal{E}_j} \left\| \mathbf{m}_j - \widetilde{\mathbf{m}}_j \right\|^2 + 2 \left(\mathbf{m}_j - \widetilde{\mathbf{m}}_j\right)^T \underbrace{\sum_{\mathbf{x}_i \in \mathcal{E}_j} \left(\mathbf{x}_i - \mathbf{m}_j\right)}_{= 0} \\ & \geq \sum_{\mathbf{x}_i \in \mathcal{E}_j} \left\|\mathbf{x}_i - \mathbf{m}_j \right\|^2 \left(= \sum_{\mathbf{x}_i \in \mathcal{E}_j} d^2(\mathbf{x}_i,\mathbf{m}_j) \right) \end{aligned}\end{split}\]

Comme \(\phi_\mathcal{E}(\mathcal{C}) \geq 0\), le processus itératif doit converger. La solution obtenue sera néanmoins un minimum local, dépendant de l’initialisation, souvent de valeur \(\phi_\mathcal{E}(\mathcal{C})\) bien plus élevée que celle correspondant au minimum global (inconnu).

La figure suivante illustre le principe de la minimisation itérative d’une fonction différentiable de deux variables. Les variables sont ici associées aux axes horizontaux et la valeur de la fonction à minimiser à l’axe vertical. A partir de l’état initial 1, un minimum (local, dans cette figure) est atteint après 3 étapes.

Minimisation itérative d’une fonction de deux variables, différentiable

Fig. 65 Minimisation itérative d’une fonction de deux variables, différentiable

Contrairement au cas illustré ci-dessus, la fonction \(\phi_\mathcal{E}(\mathcal{C})\) n’est pas différentiable (ni même continue) par rapport aux \(\mathbf{m}_j\) car un changement infinitésimal dans la position d’un centre peut provoquer un changement d’affectation de données aux centres et donc un changement significatif de la valeur de \(\phi_\mathcal{E}(\mathcal{C})\).

On appelle parfois k-moyennes (K-means) une version stochastique de la méthode des centres mobiles décrite ci-dessous (les entrées et les sorties sont les mêmes que pour l’algorithme précédent) :

  1. Initialisation aléatoire des centres \(\mathbf{m}_j\), \(1 \leq j \leq k\) ;

  2. tant que (centres non stabilisés) faire

    1. Choix aléatoire d’une des données ;

    2. Affectation de la donnée au groupe du centre le plus proche ;

    3. Recalcul des centres pour le groupe rejoint par la donnée et le groupe quitté par la donnée ;

fin tant que

Dans cette perspective, l’algorithme des centres mobiles est une variante « par lots » (batch) de K-means. Toutefois, K-means est souvent utilisé comme parfait synonyme de l’algorithme des centres mobiles ; c’est ainsi que ce terme sera employé dans la suite du cours. Vous trouverez ici un bref historique de K-means.

Initialisation de K-means : K-means++

L’initialisation des centres est une étape critique de K-means. La figure suivante illustre les résultats obtenus avec K-means à partir de trois autres initialisations aléatoires des centres :

Résultats différents avec initialisations différentes de *K-means*

Fig. 66 Résultats différents avec initialisations aléatoires différentes de K-means

On constate que ces trois solutions sont sous-optimales : pour chacune, un des groupes « naturels » (on peut en juger visuellement car les données sont bidimensionnelles et peu nombreuses, ce ne sera pas le cas en général) est divisé en deux, et deux des groupes « naturels » sont regroupés. Les valeurs de \(\phi_\mathcal{E}(\mathcal{C})\) sont plus élevées pour ces trois résultats que pour le résultat optimal illustré plus haut. Exécuter K-means plusieurs fois, à partir d’initialisations aléatoires différentes, ne garantit pas de trouver la solution optimale (ou, dans le cas général, une solution proche de la solution optimale).

Une bonne initialisation de l’algorithme K-means doit permettre d’obtenir une solution de meilleure qualité mais aussi une convergence plus rapide (c’est à dire avec moins d’itérations) vers cette solution.

Parmi les nombreux algorithmes d’initialisation de K-means nous considérerons ici K-means++ proposé dans [AV07]. L’idée de départ est assez simple : choisir les centres l’un après l’autre, suivant une loi non uniforme qui évolue après le choix de chaque centre et qui privilégie les candidats éloignés des centres déjà sélectionnés. L’algorithme d”initialisation K-means++ est :

Entrées : ensemble \(\mathcal{E}\) de \(N\) données de \(\mathbb{R}^p\) et nombre souhaité de centres \(k\) ;

Sorties : \(\mathcal{C} = \{\mathbf{c}_j, 1 \leq j \leq k\}\), à utiliser après comme centres de départ dans K-means ;

  1. \(\mathcal{C} \leftarrow\) un \(\mathbf{x}\) de \(\mathcal{E}\) choisi au hasard ;

  2. tant que (\(\|\mathcal{C}\| \le k\)) faire

    1. Sélectionner \(\mathbf{x} \in \mathcal{E}\) avec la probabilité \(\frac{d^2(\mathbf{x},\mathcal{C})}{\phi_\mathcal{E}(\mathcal{C})}\) ;

    2. \(\mathcal{C} \leftarrow \mathcal{C} \cup \{\mathbf{x}\}\) ;

fin tant que

Dans cet algorithme, nous avons noté \(d^2(\mathbf{x},\mathcal{C}) = \min_{j=1,\ldots,t} d^2(\mathbf{x},\mathbf{c}_j)\), la fonction (21) peut donc être écrite \(\phi_\mathcal{E}(\mathcal{C}) = \sum_{\mathbf{x} \in \mathcal{E}} d^2(\mathbf{x},\mathcal{C})\).

Les figures suivantes illustrent l’évolution des probabilités lors du déroulement de l’algorithme K-means++. La probabilité de sélection d’une donnée dans l’étape suivante (proportionnelle à \(d^2(\mathbf{x},\mathcal{C})\)) est d’autant plus élevée que la couleur est proche du rouge. La première coordonnée correspond à l’axe horizontal, la seconde à l’axe vertical.

Probabilités après la sélection du premier centre

Fig. 67 Probabilités après la sélection du premier centre (en bas à gauche) : \(\mathcal{C} = \left\{\left(\begin{array}[h]{c}-4,6\\8,0\end{array}\right)\right\}\)


Probabilités après la sélection du deuxième centre

Fig. 68 Probabilités après la sélection du deuxième centre (en haut à droite) : \(\mathcal{C} = \left\{\left(\begin{array}[h]{cc}-4,6 & 2,15\\8,0 & -3,45\end{array}\right)\right\}\)


Probabilités après la sélection du troisième centre

Fig. 69 Probabilités après la sélection du troisième centre (en bas à droite) : \(\mathcal{C} = \left\{\left(\begin{array}[h]{ccc}-4,6 & 2,15 & 6,32\\8,0 & -3,45 & 8,22\end{array}\right)\right\}\)


Probabilités après la sélection du quatrième centre

Fig. 70 Probabilités après la sélection du quatrième centre (en haut à gauche) : \(\mathcal{C} = \left\{\left(\begin{array}[h]{cccc}-4,6 & 2,15 & 6,32 & -8,37\\8,0 & -3,45 & 8,22 & -4,54\end{array}\right)\right\}\)

La sélection séquentielle des centres opérée par K-means++ permet de bien répartir les centres dans les données. Cela diminue le risque d’avoir, d’une part, des groupes de données sans aucun centre à proximité au démarrage de K-means et, d’autre part, plusieurs centres très proches entre eux. Ce sont ces insuffisances d’une initialisation aléatoire uniforme qui avaient provoqué la convergence vers des solutions sous-optimales dans les cas simples illustrés plus tôt.

Intérêt et limitations de K-means

Au-delà de sa simplicité conceptuelle, K-means présente d’autres intérêts qui justifient son utilisation fréquente. Il faut mentionner d’abord sa complexité algorithmique comparativement faible : \(O(t k N)\), \(t\) étant le nombre d’itérations, \(k\) le nombre de groupes (ou de centres) et \(N\) le nombre de données (observations). A première vue, la complexité serait donc linéaire en \(N\). Si on considère toutefois que le nombre \(k\) de groupes n’est pas fixé mais dépend de \(N\) (par ex. \(k \sim N^{1/3}\)) alors la complexité est légèrement sur-linéaire en \(N\) mais reste faible si on compare à d’autres méthodes. Il faut noter également que le nombre d’itérations nécessaire pour obtenir la convergence peut être fortement réduit avec une bone initialisation, utilisant par ex. K-means++. Un autre intérêt est la présence d’un unique paramètre, qui est la valeur souhaitée pour le nombre de groupes (\(k\)). En revanche, les résultats peuvent dépendre fortement de la valeur de \(k\) (surtout pour un nombre faible de groupes) et la « bonne » valeur est difficile à choisir (voir plus loin).

K-means présente également un certain nombre de limitations, dont on peut parfois s’affranchir assez facilement.

Une limitation évidente est la nécessité de disposer de données vectorielles pour pouvoir calculer les centres des groupes. Cette limitation est levée dans la méthode des k-medoids qui se satisfait de l’existence d’une métrique (au prix d’une augmentation de la complexité algorithmique).

Une limitation plus profonde est liée au critère de regroupement. Comme pour K-means chaque groupe est représenté lors du déroulement de l’algorithme par un unique prototype (son centre), la méthode n’est pas adaptée à la découverte de groupes denses et non « compacts » séparés par des régions moins denses. Des extensions qui emploient plusieurs prototypes par groupe permettent de lever dans une certaine mesure cette limitation, mais le choix le plus fréquent dans une telle situation est d’avoir recours à une méthode très différente, comme la classification spectrale ou la classification sur la base de la densité (par ex. DBSCAN et algorithmes dérivés).

La limitation de K-means aux classes de forme globalement sphérique est due à l’utilisation de la distance euclidienne usuelle et peut être facilement levée par l’emploi d’autres métriques (comme celle de Mahalanobis), au prix de l’estimation de paramètres supplémentaires (pour la distance de Mahalanobis il faut calculer une matrice de covariances empiriques par groupe).

L’utilisation de la moyenne pour calculer les centres des groupes rend K-means sensible à la présence de données atypiques ou aberrantes. Cette sensibilité peut être réduite par une estimation robuste de la moyenne, l’emploi d’un critère à minimiser plus robuste ou des extensions floues de la méthode (comme fuzzy K-means). Les extensions floues présentent par ailleurs de meilleures propriétés de convergence et permettent d’obtenir des degrés d’appartenance de chaque donnée aux différents groupes.

Enfin, si le choix du nombre \(k\) de groupes a un fort impact sur les résultats et peut rarement être fait en amont de l’application de K-means, des méthodes existent pour sélectionner une bonne valeur pour \(k\) (voir plus loin).

Méthode des K-medoids

L’emploi dans K-means des centres de gravité des groupes limite cette méthode aux données vectorielles. Dans de nombreux cas, nous souhaitons trouver des groupes dans des données pour lesquelles aucune représentation vectorielle convenable n’est disponible, par ex. des séquences de longueur variable, des sous-ensembles d’un très grand ensemble, des arbres ou des graphes. Lorsqu’une métrique (ou mesure de distance) est disponible pour ces données, il est facile d’adapter l’idée de K-means en remplaçant les centres de gravité des grupes par les « médoïdes » (medoids). Le médoïde d’un groupe est l’élément le plus « central » du groupe, c’est à dire celui pour lequel la somme des distances aux autres éléments du groupe est la plus faible.

Nous obtenons ainsi l’algorithme des K-medoids, dans lequel le seul changement par rapport à K-means est le remplacement des centres de gravité par des médoïdes :

Entrées : ensemble \(\mathcal{E}\) de \(N\) données, métrique \(d\) ;

Sorties : \(k\) groupes (clusters) disjoints \(\mathcal{E}_1,\ldots,\mathcal{E}_k\) et ensemble \(\mathcal{C}\) de leurs médoïdes ;

  1. Initialisation aléatoire : choix aléatoire de \(k\) données comme médoïdes ;

  2. tant que (medoids non stabilisés) faire

    1. Choix, pour chaque donnée, du médoïde \(\mathbf{m}_{C(l)}\) le plus proche : \(C(l) = \arg \min_{j} d(\mathbf{x}_l, \mathbf{m}_j)\) ;

    2. Constitution des groupes : \(\mathcal{E}_j\) est constitué de tous les \(\mathbf{x}_l\) qui sont plus proches de \(\mathbf{m}_j\) que de tout autre médoïde ;

    3. Recherche des médoïdes de ces (nouveaux) groupes : \(\mathbf{m}_j = \arg \min_{\mathbf{x}_l \in \mathcal{E}_j} \sum_{\mathbf{x}_p \in \mathcal{E}_j} d(\mathbf{x}_l, \mathbf{x}_p)\) ;

fin tant que

Cette version est issue de [HTF01]. Noter toutefois que l’algorithme cherche à minimiser un coût total qui est la somme des distances entre les données et les médoïdes des groupes auxquels ces données sont assignées. La première variante, Partitioning Around Medoids ([KR87]), est basée sur des essais successifs pour remplacer chaque médoïde par une autre donnée afin de réduire le coût total.

Comme pour K-means, l’initialisation a un impact important sur la qualité de la solution (du minimum local obtenu) et sur la rapidité de la convergence (le nombre d’itérations nécessaire). Il est conseillé d’employer une initialisation suivant le même principe que K-means++, c’est à dire la sélection séquentielle des centres (médoïdes ici) en maximisant l’éloignement de chaque nouveau centre des centres déjà choisis.

L’application de K-medoids à des données vectorielles, en utilisant la même métrique euclidienne que K-means, met en évidence une meilleure robustesse de K-medoids par rapport à K-means. Cette robustesse est le résultat de l’utilisation de médoïdes, moins affectés que les centres de gravité par les données excentrées.

En revanche, le remplacement des centres de gravité par des médoïdes présente un désavantage significatif en termes de complexité algorithmique. En effet, l’étape de recherche des médoïdes a pour conséquence une complexité \(O(N^2)\), alors que la complexité de K-means est plutôt \(O(N)\) (ou \(O(N^{4/3})\) si on tient compte du fait que le nombre de centres dépend de \(N\)).

Validité de la classification

Quelles que soient les données à regrouper, la plupart des algorithmes de classification automatique, dont K-means et K-medoids, convergent vers un résultat comportant en général des groupes. Quelle est la signification de ce résultat ? Y a-t-il réellement des groupes dans les données ? Leur nombre correspond ou non au nombre de groupes issu de l’algorithme ? Quel algorithme donne de meilleurs résultats sur mes données ? Nous abordons brièvement ces questions dans la suite.

Choix du nombre de groupes \(k\)

Les méthodes de partitionnement que nous avons vues jusqu’ici ne modifient pas le nombre de groupes (\(k\)) fourni en entrée de l’algorithme. Comment choisir alors le « bon » nombre de groupes, vu que des informations permettant de le définir a priori sont très rarement disponibles ?

Nous avons également constaté que la valeur minimale que pouvait atteindre le critère (21) diminuait avec l’augmentation de \(k\) (et peut même arriver à \(0\) pour \(k = N\)), il n’est donc pas envisageable d’utiliser ce critère pour choisir entre plusieurs valeurs de \(k\). Parmi les différentes approches employées pour le choix de \(k\) nous pouvons mentionner :

  1. La méthode du « coude ». Il est nécessaire d’appliquer l’algorithme de classification automatique (par ex. K-means avec initialisation par K-means++ ) pour un ensemble de valeurs du nombre de groupes \(k\) et de représenter sur un graphique les valeurs minimales atteintes par \(\phi_\mathcal{E}(\mathcal{C})\) (voir Fig. 71). La valeur retenue pour \(k\) est celle qui marque le début d’un pallier : pour des valeurs inférieures la qualité de regroupement est nettement moins bonne, alors que pour des valeurs supérieures la qualité ne s’améliore pas sensiblement. Noter toutefois que si la plage de valeurs de \(k\) explorées est suffisamment large il est fréquent d’y trouver plusieurs « coudes ».

    Choix du nombre de groupes avec la méthode du coude

    Fig. 71 Choix du nombre \(k\) de groupes avec la méthode du coude

  2. La reformulation du problème de classification automatique dans un cadre probabiliste. Cela permet de choisir un critère mieux fondé, par exemple un critère d’information comme AIC (Akaike) ou BIC (Bayes) pour identifier la valeur optimale pour \(k\), voir le chapitre sur l’estimation de densité.

  3. L’étude de la stabilité des résultats. Des résultats relativement récents (voir par ex. [ST10]) proposent la stabilité comme critère de choix de modèle : par ex., une valeur de \(k\) est meilleure si les groupes obtenus sont plus stables à l’initialisation aléatoire de l’algorithme.

Les méthodes mentionnées permettent de choisir a posteriori une bonne valeur pour le nombre \(k\) de groupes au prix de calculs plus lourds, car l’algorithme de classification doit être appliqué une fois (ou même plusieurs fois afin d’évaluer la stabilité) pour chaque valeur candidate.

Mesures de cohérence entre classifications

Les résultats obtenus par classification automatique sont d’autant plus significatifs que deux méthodes de classification différentes, ou deux initialisations différentes pour une même méthode sensible à l’initialisation, produisent (approximativement) les mêmes groupes. Pour mesurer la cohérence entre les résultats de deux applications différentes d’algorithmes de classification il est possible d’employer différentes mesures, par ex. (voir [WW07] [VEB10]) :

  1. L’indice de Rand ajusté. Une classification \(\mathcal{C}\) est définie comme l’ensemble des paires d’observations qui sont dans un même groupe, parmi toutes les paires possibles. On note par \(n_{11}\) le nombre de paires qui sont dans un même groupe suivant \(\mathcal{C}\) et \(\mathcal{C}'\), par \(n_{00}\) le nombre de paires qui sont dans des groupes différents suivant \(\mathcal{C}\) et \(\mathcal{C}'\), par \(n_{10}\) le nombre de paires dans un même groupe suivant \(\mathcal{C}\) et groupes différents suivant \(\mathcal{C}'\), et enfin par \(n_{01}\) le nombre de paires dans un même groupe suivant \(\mathcal{C}'\) et groupes différents suivant \(\mathcal{C}\).

    Si \(\mathcal{C}, \mathcal{C}'\) sont deux classifications différentes, l’indice de Rand est la probabilité pour que les deux classifications soient en accord pour une paire de données choisie au hasard. Cela peut s’exprimer sous la forme suivante

    (22)\[\mathcal{R}(\mathcal{C},\mathcal{C}') = \frac{2 (n_{11} + n_{00})}{N (N - 1)}, \ \ 0 \leq \mathcal{R} \leq 1\]

    L’indice de Rand ajusté est défini par

    (23)\[\mathcal{R}_{adj}(\mathcal{C},\mathcal{C}') = \frac{\mathcal{R} - E(\mathcal{R})}{\max(\mathcal{R}) - E(\mathcal{R})}\]

    \(E(\mathcal{R})\) est l’espérance de l’indice de Rand et \(\max(\mathcal{R})\) la valeur maximale qu’il peut atteindre. Nous nous servirons de l’indice de Rand ajusté dans la séance de travaux pratiques.

  2. L’indice de Jaccard, défini par
    \[\mathcal{I}(\mathcal{C},\mathcal{C}') = \frac{n_{11}}{n_{11}+n_{10}+n_{01}}\]
  3. L’information mutuelle normalisée, pour laquelle plusieurs variantes existent, par ex.

    \[NMI(\mathcal{C},\mathcal{C}') = \frac{I(\mathcal{C},\mathcal{C}')}{\sqrt{H(\mathcal{C}) H(\mathcal{C}')}}\]

    \(I(\mathcal{C},\mathcal{C}')\) est l’information mutuelle entre \(\mathcal{C},\mathcal{C}'\) et \(H(\mathcal{C})\), \(H(\mathcal{C}')\) les entropies respectives des deux partitionnements.

L’indice de Rand ajusté et l’information mutuelle normalisée sont disponibles, avec d’autres mesures, dans Scikit-learn.

Validation interne

La qualité des résultats obtenus par une méthode de classification automatique peut être évaluée en examinant dans quelle mesure les groupes identifiés sont bien « ajustés » aux données. Par exemple, dans quelle mesure les données d’une même groupe sont plus proches entre elles que des données d’autres groupes.

De nombreux critères ont été proposés pour cela, comme le coefficient de silhouette (mesure à quel point les groupes sont « compacts » et bien séparés entre eux), la statistique modifiée de Hubert (mesure l’alignement entre distance et partitionnement), l’indice Davies-Bouldin (qui peut être vu comme un rapport entre inerties intra-groupes et inerties inter-groupes), etc. Plusieurs critères de ce type sont disponibles dans Scikit-learn. En général, ces critères ne peuvent pas être optimisés directement par un algorithme de classification automatique.

Il est important de noter que ces critères reflètent à la fois les propriétés des données elles-mêmes (forment-elles naturellement des groupes compacts et bien séparés) et l’ajustement aux données de la solution trouvée par l’algorithme de classification automatique. Cette caractéristique doit être prise en compte lorsqu’on compare les valeurs d’un même critère sur deux problèmes différents (c’est à dire avec des données différentes).

Nous pouvons remarquer que les critères (ou fonctions de coût) minimisés par les différentes méthodes de classification automatique (comme (21) pour K-means) reflètent également l’ajustement des groupes aux données mais ne permettent pas de comparer entre eux les résultats sur deux problèmes différents, ni de comparer à une référence.

Validation relative

Pour obtenir les meilleurs résultats de classification automatique sur un ensemble de données il est nécessaire d’appliquer plusieurs méthodes de classification différentes et de choisir ensuite celle qui mène aux meilleurs résultats. Les critères d’évaluation mentionnés dans la section précédente sont bien adaptés à ce type de comparaison car sur un même ensemble de données les différences mesurées concernent exclusivement l’ajustement des groupes trouvés aux données. En revanche, ces critères peuvent être plus ou moins adaptés aux données à traiter, par ex. le coefficient de silhouette n’est pas le critère adéquat pour des données qui se regroupent sur la base de la densité et non de la compacité.

Validation externe

Dans certains cas il est possible que des connaissances a priori concernant le regroupement ou non d’une partie des données soient disponibles. Par exemple, le fait que certaines données appartiennent à certains groupes (information non disponible pour tous les groupes et encore moins pour toutes les données), ou que certaines données doivent appartenir à un même groupe (contraintes must-link) ou à des groupes différents (contraintes cannot-link). Ces connaissances ne sont pas nécessairement directement exploitables dans un algorithme de classification (par ex. dans une fonctionnelle à minimiser), bien que des avancées récentes ont été faites dans le domaine de la classification semi-supervisée (semi-supervised clustering, voir par ex. [Bai13]). En revanche, ces connaissances peuvent toujours être employées à l’issue de la classification automatique pour évaluer les résultats obtenus.

Classification ascendante hiérarchique

Parmi les méthodes qui visent à obtenir une hiérarchie de regroupements, nous examinons brièvement dans la suite la classification ascendante hiérarchique (CAH). Par rapport à un simple partitionnement des données, une hiérarchie de regroupements fournit une information plus riche concernant la structure de similarité des données. En effet, elle permet d’examiner l’ordre des agrégations de groupes, les rapports des similarités entre groupes, etc.

La classification ascendante procède par agrégations successives de groupes. A partir de la hiérarchie de groupes résultante il est possible d’observer l’ordre des agrégations de groupes, d’examiner les rapports des similarités entre groupes, ainsi que d’obtenir plusieurs partitionnements à des niveaux de granularité différents.

L’exemple suivant montre la hiérarchie de groupes (« dendrogramme ») obtenue par agrégations successives à partir d’un petit ensemble de données bidimensionnelles :

Dendrogramme obtenu à partir de données 2D

Fig. 72 Dendrogramme (à droite) obtenu à partir des données 2D (à gauche)

Le principe de la classification ascendante est de regrouper (ou d’agréger), à chaque itération, les données et/ou les groupes les plus proches qui n’ont pas encore été regroupé(e)s.

Considérons un ensemble de données \(\mathcal{E} \subset \mathcal{X}\) et une distance \(d_\mathcal{X} : \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}^+\). Les données individuelles peuvent être comparées grâce à cette distance. Mais, une fois deux données agrégées dans un groupe (noté \(h_p\), par ex.), comment mesurer la similarité de ce groupe à une autre donnée ou à un autre groupe (\(h_q\), par ex.) ? Il est nécessaire de définir une nouvelle mesure pour la dissimilarité entre groupes (ou entre une donnée et un groupe). Cette mesure est en général appelée « indice d’agrégation » et plusieurs choix fréquents sont :

  1. Le lien minimum (single linkage) : l’indice d’agrégation entre deux groupes \(h_p,h_q\) est la valeur la plus faible des distances entre une donnée du premier groupe et une donnée du second. Formellement, \(\delta_s(h_p,h_q) = \min_{x_i \in h_p, x_j \in h_q} d_\mathcal{X}(x_i,x_j)\). Il suffit donc que deux groupes contiennent deux éléments proches pour que la valeur de l’indice soit faible :

Illustration du lien minimum

Fig. 73 Illustration du lien minimum (single linkage)

  1. Le lien maximum (complete linkage) : l’indice d’agrégation entre deux groupes \(h_p,h_q\) est la valeur la plus élevée des distances entre une donnée du premier groupe et une donnée du second (appelée parfois « diamètre » de l’agrégat). Formellement, \(\delta_s(h_p,h_q) = \max_{x_i \in h_p, x_j \in h_q} d_\mathcal{X}(x_i,x_j)\). Il est donc nécessaire que tous les éléments d’un des groupes soient proches de tous les éléments de l’autre groupe pour que la valeur de l’indice soit faible :

Illustration du lien maximum

Fig. 74 Illustration de l’agrégation utilisant le lien maximum (complete linkage)

  1. Le lien moyen (average linkage) : l’indice d’agrégation entre deux groupes \(h_p,h_q\) est la moyenne des distances entre une donnée du premier groupe et une donnée du second. Formellement, \(\delta_s(h_p,h_q) = \frac{1}{\|h_p\|\cdot\|h_q\|} \sum_{x_i \in h_p, x_j \in h_q} d_\mathcal{X}(x_i,x_j)\) :

Illustration du lien moyen

Fig. 75 Illustration de l’agrégation utilisant le lien moyen (complete linkage)

  1. Les illustrations précédentes représentent des données de \(\mathbb{R}^2\), mais les définitions des trois indices montrent bien que la seule exigence est l’existence d’une distance entre données individuelles et non d’une structure d’espace vectoriel.

    L’indice de Ward, en revanche, est défini uniquement pour des données vectorielles. La valeur de l’indice correspond à l’augmentation de l’inertie intra-groupe lors de l’agrégation de deux groupes. Elle est la différence entre l’inertie intra-groupe de l’agrégat et la somme des inerties intra-groupe des deux groupes agrégés. Formellement, \(\delta_s(h_p,h_q) = \frac{\|h_p\|\cdot\|h_q\|}{\|h_p\|+\|h_q\|} d^2_\mathcal{X}(\mathbf{m}_p,\mathbf{m}_q)\).

Illustration de l'indice de Ward

Fig. 76 Illustration de l’agrégation utilisant l’indice de Ward

Contrairement à l’algorithme K-means, qui exige des données vectorielles, la CAH se satisfait, pour de nombreux indices d’agrégation (par ex. lien minimum, lien maximum, lien moyen) d’une simple distance définie sur l’espace des données. Il est donc possible de l’employer directement pour des données qui n’ont pas de représentation vectorielle naturelle.

L’algorithme de CAH est très simple : on agrège, lors d’itérations successives, les groupes entre lesquels l’indice d’agrégation choisi a la valeur minimale.

Entrées : \(\mathcal{E}\) de \(N\) données de \(\mathcal{X}\) muni de la distance \(d_\mathcal{X}\), indice d’agrégation \(\delta_s\) ;

Sorties : Hiérarchie de groupes (dendrogramme) ;

  1. Chaque donnée définit un groupe ;

  2. tant que (nombre de groupes \(> 1\)) faire

    1. Calcul des valeurs de l’indice d’agrégation entre tous les groupes issus de l’itération précédente ;

    2. Regroupement des 2 groupes ayant la plus petite valeur de l’indice d’agrégation ;

fin tant que

Dans la hiérarchie de regroupements (dendrogramme) illustrée plus haut, la hauteur de chaque palier (qui représente une agrégation) est égale à la valeur de l’indice d’agrégation entre les groupes agrégés.

La complexité algorithmique pour un choix non contraint de l’indice d’agrégation est \(O(N^2 \log N)\). En revanche, pour l’indice du lien minimum et l’indice du lien maximum existent des algorithmes de complexité \(O(N^2)\), complexité qui reste néanmoins élevée. Une pratique assez fréquente est d’appliquer d’abord l’algorithme K-means avec une valeur élevée de \(k\) (mais néanmoins \(k \ll N\)), ensuite d’utiliser la classification ascendante hiérarchique pour regrouper les \(k\) groupes issus de K-means.

Le choix de l’indice d’agrégation a un impact parfois significatif sur les résultats obtenus. Au-delà des contraintes liées à la nature des données (éléments d’un espace vectoriel ou simplement métrique), une bonne compréhension du problème doit aider ce choix. Par exemple, l’utilisation du lien minimum mène à un résultat qui se rapproche du regroupement par densité, alors que l’emploi du lien maximum (ou de l’indice de Ward) produit plutôt des groupes « compacts » séparés entre eux. En revanche, le lien minimum mène souvent à des arbres en escalier, déséquilibrés et peu exploitables, comme dans l’exemple suivant :

Illustration de l'effet escalier

Fig. 77 Illustration de l’effet escalier obtenu avec l’indice du lien minimum

Quelle que soit la hauteur à laquelle on coupe un tel arbre, le partitionnement obtenu n’est pas exploitable car le résultat correspond à un grand groupe plus un ensemble de données individuelles.

Classification spectrale

« Si vous pensez avez compris la classification spectrale c’est qu’on ne vous l’a pas assez bien expliquée. » (Un auditeur)

La classification spectrale vise à obtenir un partitionnement des données (ou observations) en groupes à partir d’une représentation de ces données sous la forme d’un graphe de similarités. Dans un tel graphe, chaque sommet correspond à une donnée (ou observation) et chaque arête qui relie deux observations est pondérée par la similarité entre ces observations.

Plus précisément, soit un ensemble \(\mathcal{E} = \{\mathbf{x}_i, 1 \leq i \leq N\}\) de \(N\) données décrites par une matrice de similarités \(\mathbf{S}\) telle que l’élément \(s_{ij} \geq 0\) soit la similarité entre \(\mathbf{x}_i\) et \(\mathbf{x}_j\).

L’objectif de la classification automatique est de répartir les \(N\) données de \(\mathcal{E}\) en \(k\) groupes disjoints \(\mathcal{E}_1,\ldots,\mathcal{E}_k\) (inconnus a priori) tels que les similarités soient fortes intra-groupe et faibles entre les groupes.

La méthode de classification spectrale procède en plusieurs étapes :

  1. A partir de \(\mathbf{S}\) est obtenu un graphe de similarités. Plusieurs méthodes alternatives de construction d’un tel graphe sont examinées dans la suite.

  2. A l’aide du graphe, une représentation vectorielle « convenable » (du point de vue de la classification automatique) des données est obtenue.

  3. Un algorithme de classification de base (comme K-means) est enfin appliqué pour obtenir les \(k\) groupes disjoints \(\mathcal{E}_1,\ldots,\mathcal{E}_k\).

Afin d’illustrer le fonctionnement de la classification spectrale, nous nous servirons dans la suite de l’exemple donné dans [Lux07]. Les données sont des vecteurs dans l’espace bidimensionnel, illustrés dans la figure suivante par des étoiles bleues.

Données pour l'exemple de classification spectrale

Fig. 78 Données pour l’exemple de classification spectrale (issu de [Lux07])

A partir de ces données vectorielles, une matrice de similarités est obtenue et ensuite un graphe de similarités est déduit. Nous examinons plus loin les méthodes permettant de calculer des similarités à partir des distances et d’obtenir un graphe à partir de la matrice de similarités. Dans l’exemple considéré nous pouvons illustrer un graphe de similarités par la figure suivante.

Graphe de similarités pour l'exemple

Fig. 79 Graphe de similarités pour l’exemple (issu de [Lux07])

Avant d’aller plus loin, il est nécessaire d’introduire (ou rappeler) quelques notions concernant les graphes.

Un graphe non orienté est un couple \(G = (V,E)\), où \(V = \{v_1, \ldots, v_N\}\) est l’ensemble des (\(N\)) sommets et \(E\) l’ensemble des arêtes. Voir l’exemple de graphe de la figure précédente, où les sommets sont les données (vecteurs de \(\mathbb{R}^2\) dans ce cas).

Les arêtes peuvent être pondérées : soit \(w_{ij} \geq 0\) le poids de l’arête qui lie les sommets \(v_i\) et \(v_j\), alors \(w_{ij} = 0\) si et seulement si (ssi) il n’y a pas d’arête entre les sommets \(v_i\) et \(v_j\). La matrice des poids est \(\mathbf{W}\) d’élément générique \(w_{ij}\).

Le degré du sommet \(v_i\) est la somme des poids des arêtes qui relient ce sommet à d’autres, \(d_i = \sum_{j=1}^N w_{ij}\). On peut considérer une matrice diagonale des degrés \(\mathbf{D}\) qui a \(\{d_1, \ldots, d_N\}\) sur la diagonale principale.

Une chaîne ou un chemin \(\mu(v_i,v_j)\) entre les sommets \(v_i\) et \(v_j\) est une suite d’arêtes de \(E\) permettant de relier les deux sommets.

Un sous-ensemble \(C \subset V\) forme une composante connexe du graphe \(G = (V,E)\) si

  • quels que soient deux sommets de \(C\), il existe une chaîne de \(G\) les reliant (\(\forall v_i,v_j \in C\), \(\exists \mu(v_i,v_j)\)) et

  • pour tout sommet \(v_i\) de \(C\) et toute chaîne avec des arêtes de \(E\) et partant de \(v_i\), le sommet terminal de la chaîne est également dans \(C\) (\(\forall v_i \in C, \forall \mu(v_i,v_p)\), \(v_p \in C\)).

Le vecteur indicateur d’un sous-ensemble \(C\) de \(V\) est le vecteur \(\mathbf{c}\) à \(N\) composantes binaires tel que \(c_i = 1\) si le sommet \(v_i \in C\) et \(c_i = 0\) sinon. La même définition s’applique lorsque \(C\) est une composante connexe.

Pour un graphe non orienté à arêtes pondérées, la matrice laplacienne (non normalisée) est définie par \(\mathbf{L} = \mathbf{D} - \mathbf{W}\). La matrice laplacienne possède plusieurs propriétés importantes :

  • pour tout vecteur \(\mathbf{x} \in \mathbb{R}^N\), nous avons \(\mathbf{x}^T \mathbf{L} \mathbf{x} = \frac{1}{2} \sum_{i,j=1}^N w_{ij} (x_i - x_j)^2\) ;

  • la matrice laplacienne \(\mathbf{L}\) est symétrique (\(\mathbf{L} = \mathbf{L}^T\)) et semi-définie positive (\(\forall \mathbf{x} \in \mathbb{R}^N\), \(\mathbf{x}^T \mathbf{L} \mathbf{x} \geq 0\));

  • la matrice laplacienne \(\mathbf{L}\) a \(N\) valeurs propres réelles \(0 = \lambda_1 \leq \lambda_2 \leq \ldots \leq \lambda_N\) ;

  • la valeur propre la plus faible de \(\mathbf{L}\) est 0, le vecteur propre correspondant étant le vecteur constant \(\mathbf{1}\) (toutes ses composantes sont égales à \(1\)) ; cela s’explique facilement par la définition même des degrés des sommets ;

  • la multiplicité de la valeur propre 0 de \(\mathbf{L}\) est égale au nombre de composantes connexes de \(G\) ;

  • les vecteurs indicateurs de ces composantes connexes sont des vecteurs propres de \(\mathbf{L}\) pour la valeur propre 0.

Plusieurs méthodes peuvent être employées pour normaliser la matrice laplacienne d’un graphe, notamment :

  1. La normalisation symétrique : \(\mathbf{L}_{sym} = \mathbf{D}^{-1/2} \mathbf{L} \mathbf{D}^{-1/2} \left(= \mathbf{I}_N - \mathbf{D}^{-1/2} \mathbf{W} \mathbf{D}^{-1/2}\right)\).

  2. La normalisation « marche aléatoire » (random walk), où le résultat n’est plus une matrice symétrique : \(\mathbf{L}_{rw} = \mathbf{D}^{-1} \mathbf{L} \left(= \mathbf{I}_N - \mathbf{D}^{-1} \mathbf{W}\right)\).

Algorithme de classification spectrale

Nous pouvons maintenant détailler les étapes 2 et 3 de l’algorithme de classification spectrale :

  1. A partir de \(\mathbf{S}\) est produit un graphe de similarités \(G\), avec les pondérations \(w_{ij} = s_{ij}\).

  2. A l’aide du graphe, une représentation vectorielle « convenable » (du point de vue de la classification automatique) des données est obtenue par la méthode suivante :

  1. Calculer la matrice laplacienne \(\mathbf{L}\) ;

  2. Calculer les \(k\) vecteurs propres \(\mathbf{u}_1, \ldots, \mathbf{u}_k\) associés aux \(k\) plus petites valeurs propres de \(\mathbf{L} \mathbf{u} = \lambda \mathbf{D} \mathbf{u}\) ;

  3. Soit \(\mathbf{U}_k\) la matrice dont les colonnes sont les vecteurs \(\mathbf{u}_1, \ldots, \mathbf{u}_k\) ;

  4. Soit \(\mathbf{y}_i \in \mathbb{R}^k, 1 \leq i \leq N\), les \(N\) lignes de la matrice \(\mathbf{U}_k\), chaque donnée (observation) \(i\) est représentée par le vecteur \(\mathbf{y}_i\).

  1. Un algorithme de classification de base (comme K-means) est enfin appliqué pour obtenir les \(k\) groupes disjoints \(\mathcal{E}_1,\ldots,\mathcal{E}_k\) :

  1. Appliquer K-means aux \(N\) vecteurs \(\mathbf{y}_i\) pour obtenir les groupes \(\mathcal{C}_1, \ldots, \mathcal{C}_k\) ;

  2. Pour \(1 \leq j \leq k\), \(\mathcal{E}_j = \{\mathbf{x}_i | \mathbf{y}_i \in \mathcal{C}_j\}\).

Cet algorithme peut être justifié de plusieurs façons différentes. Nous mentionnons dans la suite deux des justifications les plus courantes, sans toutefois détailler complètement les formulations mathématiques correspondantes (qui peuvent être consultées par ex. dans [Lux07] et ses références).

Coupe normalisée d’un graphe

Considérons un graphe non orienté à arêtes pondérées \(G = (V, E)\), que nous souhaitons découper en \(k\) groupes (sous-ensembles) disjoints \(V_1, \ldots, V_k\) de sommets.

La coupe minimum est obtenue en choisissant les \(k\) groupes de sommets de façon à minimiser la somme des poids des arêtes qui relient des sommets appartenant à des groupes différents. Plus précisément, est minimisée la quantité \(\sum_{s=1}^k \omega(V_s, V-V_s)\), où \(\omega(V_s, V_t) = \sum_{v_i \in V_s, v_j \in V_t} w_{ij}\) et \(V-V_s\) est l’ensemble de sommets de \(V\) qui ne sont pas dans \(V_s\). Malheureusement, bien souvent la coupe minimum sépare le graphe en un nombre limité \(l \ll k\) de grands sous-graphes plus quelques sommets isolés, solution jugée en général peu satisfaisante.

Avec la coupe normalisée (normalized cut, Ncut), les sous-ensembles de sommets sont plus équilibrés. La coupe normalisée est obtenue en minimisant \(\sum_{s=1}^k \frac{\omega(V_s, V-V_s)}{\textrm{vol}(V_s)}\), où le « volume » d’une composante \(V_s\) est la somme des degrés des sommets de cette composante, \(\textrm{vol}(V_s) = \sum_{v_i \in V_s} d_i\). Malheureusement, si la solution obtenue est en général plus satisfaisante, ce problème de minimisation est NP-difficile !

La classification spectrale reformule le problème de minimisation en utilisant un vecteur indicateur pour chaque sous-ensemble de sommets, avec une « relaxation » continue du problème discret : les composantes du vecteur indicateur ne sont pas binaires mais ont des valeurs dans l’intervalle \([0, 1]\). Considérons la matrice \(\mathbf{U}_k\) dont les colonnes sont les vecteurs caractéristiques \(\mathbf{u}_1, \ldots, \mathbf{u}_k\) des \(k\) groupes recherchés et chacune des \(N\) lignes correspond à une des observations.

Avec des groupes disjoints, chacune des observations appartient à un seul groupe et donc le vecteur \(\mathbf{y}_i\) qui représente l’observation \(i\) a une seule composante égale à 1 (la composante \(y_{ij}\) si l’observation \(i\) fait partie du groupe \(j\)), alors que les autres sont égales à 0. Avec la relaxation continue, ces composantes du vecteur ont toutes des valeurs dans l’intervalle \([0, 1]\) ; on peut les voir comme des degrés d’appartenance de cette observation aux différents groupes.

Avec cette relaxation continue, la coupe normalisée est ainsi traduite en un problème de minimisation sous contraintes d’égalité. Les contraintes d’égalité consistent à imposer, pour toute ligne de \(\mathbf{U}_k\), que la somme des composantes soit égale à 1 : \(\forall i \in \{1, \ldots, N\}, \sum_{j=1}^k y_{ij} = 1\).

Pour résoudre ce nouveau problème de minimisation, il faut trouver les \(k\) plus petites valeurs propres de \(\mathbf{L} \mathbf{u} = \lambda \mathbf{D} \mathbf{u}\) et les vecteurs propres associés (qui seront les colonnes de \(\mathbf{U}_k\)). En l’absence de sommets isolés, \(\mathbf{L} \mathbf{u} = \lambda \mathbf{D} \mathbf{u}\) devient \(\mathbf{D}^{-1} \mathbf{L} \mathbf{u} = \lambda \mathbf{u}\), ou encore \(\mathbf{L}_{rw} \mathbf{u} = \lambda \mathbf{u}\).

Pour revenir à des vecteurs caractéristiques binaires, dans la dernière étape de la classification spectrale un algorithme de classification automatique de base (comme K-means) est appliqué aux lignes de \(\mathbf{U}_k\) et ainsi chaque donnée est affectée à un groupe (chaque sommet est affecté à une composante du graphe).

Marche aléatoire sur un graphe

On arrive au même problème de coupe normalisée en cherchant à identifier des composantes du graphe telles qu’une marche aléatoire a peu d’opportunités de passer d’une composante à une autre. Une marche aléatoire sur un graphe est une succession de sauts aléatoires d’un sommet du graphe à un autre sommet, suivant des probabilités de transition. La probabilité de passer du sommet \(i\) au sommet \(j\) (probabilité de transition \(i \rightarrow j\)) est le rapport entre le poids de l’arête (ici, similarité entre les sommets \(i\) et \(j\)) et la somme des poids des arêtes partant du sommet \(i\) (le degré de ce sommet) : \(p_{ij} = w_{ij} / d_i\).

Avec ces probabilités de transition, la matrice de transition est \(\mathbf{P} = \mathbf{D}^{-1} \mathbf{W}\), c’est à dire \(\mathbf{P} = \mathbf{I}_N - \mathbf{L}_{rw}\). Alors, \(\lambda\) est une valeur propre de \(\mathbf{L}_{rw}\) avec comme vecteur propre \(\mathbf{u}\) si et seulement si \(1-\lambda\) est une valeur propre de \(\mathbf{P}\) avec le même vecteur propre \(\mathbf{u}\).

On peut montrer que trouver \(V_s \subset V\) pour minimiser la probabilité de passage de \(V_s\) à \(V - V_s\) et de \(V - V_s\) à \(V_s\) est équivalent au problème de la coupe normalisée. La classification spectrale fournit une solution en général approximative au problème de la coupe normalisée et donc également à l’identification de composantes d’un graphe telles qu’une marche aléatoire a peu d’opportunités de passer d’une composante à une autre.

Deux exemples numériques simples

Deux exemples numériques liés entre eux peuvent aider à mieux comprendre comment est obtenue la représentation vectorielle « convenable » (du point de vue de la classification automatique) des données. Considérons les 5 données (observations) bidimensionnelles qui sont les lignes de la matrice suivante :

\[\begin{split}\left[\begin{array}{rr} 0 & 0 \\ 1 & 0 \\ 2 & 0 \\ 2 & 3 \\ 0 & 3 \end{array}\right]\end{split}\]

Ces données sont représentées aussi dans la figure ci-dessous :

Données pour l'exemple numérique 1

Fig. 80 Données pour l’exemple numérique 1 de classification spectrale

En faisant attention au fait que l’échelle de l’abscisse n’est pas tout à fait la même que celle de l’ordonnée, nous pouvons constater visuellement que le choix de 2 groupes ({A, B, C} et {D, E}) serait naturel dans ce cas, suivi par le choix de 3 groupes ({A, B, C}, {D}, {E}).

Le deuxième exemple est issu du premier en écartant beaucoup plus les observations D et E de A, B et C, au point où cela peut produire la déconnexion entre deux composantes du graphe :

\[\begin{split}\left[\begin{array}{rr} 0 & 0 \\ 1 & 0 \\ 2 & 0 \\ 2 & \mathbf{10} \\ 0 & \mathbf{10} \end{array}\right]\end{split}\]

La figure suivante illustre ce deuxième exemple (attention, l’échelle verticale est très différente) :

Données pour l'exemple numérique 2

Fig. 81 Données pour l’exemple numérique 2 de classification spectrale

Dans ce cas, par la visualisation on constate que le choix naturel est celui de 2 groupes : {A, B, C} et {D, E}.

Nous utilisons la distance euclidienne classique pour mesurer la distance entre les observations. La matrice des similarités est obtenue à partir du noyau gaussien de variance 1 appliqué à la distance euclidienne. Les matrices de similarités résultantes (valeurs arrondies à 6 décimales) sont alors :

  • Exemple 1 :

    \[\begin{split}\left[\begin{array}{rrrrr} 1.000000 & 0.367879 & 0.018316 & 0.000002 & 0.000123 \\ 0.367879 & 1.000000 & 0.367879 & 0.000045 & 0.000045 \\ 0.018316 & 0.367879 & 1.000000 & 0.000123 & 0.000002 \\ 0.000002 & 0.000045 & 0.000123 & 1.000000 & 0.018316 \\ 0.000123 & 0.000045 & 0.000002 & 0.018316 & 1.000000 \end{array}\right]\end{split}\]
  • Exemple 2 :

    \[\begin{split}\left[\begin{array}{rrrrr} 1.000000 & 0.367879 & 0.018316 & \mathbf{0.000000} & \mathbf{0.000000} \\ 0.367879 & 1.000000 & 0.367879 & \mathbf{0.000000} & \mathbf{0.000000} \\ 0.018316 & 0.367879 & 1.000000 & \mathbf{0.000000} & \mathbf{0.000000} \\ 0.\mathbf{000000} & \mathbf{0.000000} & \mathbf{0.000000} & 1.000000 & 0.018316 \\ 0.\mathbf{000000} & \mathbf{0.000000} & \mathbf{0.000000} & 0.018316 & 1.000000 \end{array}\right]\end{split}\]

On observe que le graphe associé au premier exemple peut être considéré complètement connecté car toutes les similarités sont strictement supérieures à 0. En revanche, pour le second exemple, avec la précision de représentation considérée, les similarités entre les observations de l’ensemble {A, B, C} et celles de l’ensemble {D, E} sont nulles (valeurs représentées en gras dans la matrice), le graphe associé possède donc deux composantes connexes, {A, B, C} et {D, E}.

Les matrices laplaciennes normalisées asymétriques \(\mathbf{L}_{rw} = \mathbf{D}^{-1} \mathbf{L} \left(= \mathbf{I}_N - \mathbf{D}^{-1} \mathbf{W}\right)\) résultantes (valeurs arrondies à 6 décimales) sont :

  • Exemple 1 :

    \[\begin{split}\left[\begin{array}{rrrrr} 1.000000 & -0.952264 & -0.047410 & -0.000006 & -0.000319 \\ -0.499938 & 1.000000 & -0.499938 & -0.000062 & -0.000062 \\ -0.047410 & -0.952264 & 1.000000 & -0.000319 & -0.000006 \\ -0.000122 & -0.002456 & -0.006676 & 1.000000 & -0.990746 \\ -0.006676 & -0.002456 & -0.000122 & -0.990746 & 1.000000 \end{array}\right]\end{split}\]
  • Exemple 2 :

    \[\begin{split}\left[\begin{array}{rrrrr} 1.000000 & -0.952574 & -0.047426 & \mathbf{0.000000} & \mathbf{0.000000} \\ -0.500000 & 1.000000 & -0.500000 & \mathbf{0.000000} & \mathbf{0.000000} \\ -0.047426 & -0.952574 & 1.000000 & \mathbf{0.000000} & \mathbf{0.000000} \\ \mathbf{0.000000} & \mathbf{0.000000} & \mathbf{0.000000} & 1.000000 & -1.000000 \\ \mathbf{0.000000} & \mathbf{0.000000} & \mathbf{0.000000} & -1.000000 & 1.000000 \end{array}\right]\end{split}\]

Pour le premier exemple, les valeurs propres de \(\mathbf{L}_{rw}\) en ordre croissant sont 0.0, 0.0094, 1.0474, 1.9523 et 1.9907. Pour obtenir une classification en 2 groupes nous nous intéressons aux deux plus petites valeurs propres de \(\mathbf{L}_{rw}\), qui sont ici 0.0 et 0.0094. Les colonnes de la matrice \(\mathbf{U}_2\) suivante sont les 2 vecteurs propres correspondants :

\[\begin{split}\left[\begin{array}{rr} 0.447214 & -0.017287 \\ 0.447214 & -0.017287 \\ 0.447214 & -0.017287 \\ 0.447214 & 0.706789 \\ 0.447214 & 0.706789 \end{array}\right]\end{split}\]

Les nouvelles représentations des données pour l’exemple 1 sont alors les lignes de la matrice \(\mathbf{U}_2\). L’illustration suivante montre que les projections sur ces 2 vecteurs propres des observations A, B et C sont confondues (au niveau de précision numérique choisi) et différentes des projections de D et E qui sont également confondues :

Nouvelles représentations des données pour l'exemple numérique 1

Fig. 82 Nouvelles représentations des données pour l’exemple numérique 1

Nous remarquerons que le vecteur propre associé à la valeur propre 0 a toutes les composantes identiques et ne présente donc aucune utilité pour distinguer entre les (projections des) observations.

Le choix de trois groupes aurait été bien moins naturel car la troisième valeur propre la plus petite est 1.0474, nettement supérieure à la deuxième qui est 0.0094.

Pour le deuxième exemple, les valeurs propres de \(\mathbf{L}_{rw}\) en ordre croissant sont 0.0, 0.0, 1.0474, 1.9525 et 2.0000. La valeur propre 0 est donc de multiplicité 2, ce qui indique qu’il y a bien 2 composantes connexes dans le graphe. La matrice \(\mathbf{U}_2\) dont les colonnes sont les 2 vecteurs propres correspondant à la même valeur propre 0 (de multiplicité 2) est :

\[\begin{split}\left[\begin{array}{rr} 0.577350 & 0.000000 \\ 0.577350 & 0.000000 \\ 0.577350 & 0.000000 \\ 0.000000 & 0.707107 \\ 0.000000 & 0.707107 \end{array}\right]\end{split}\]

Il faut noter que le choix exact des deux vecteurs propres est arbitraire dans le sous-espace propre à 2 dimensions correspondant à la valeur propre 0 (la dimension de l’espace propre est égale à la multiplicité de la valeur propre). Les seules contraintes sont l’orthogonalité de ces deux vecteurs propres et le fait que chacun a une norme égale à 1.

Les nouvelles représentations des données pour l’exemple 2 sont alors les lignes de la matrice \(\mathbf{U}_2\). L’illustration suivante montre que pour ce deuxième exemple aussi les projections sur ces 2 vecteurs propres des observations A, B et C sont confondues et différentes des projections de D et E qui sont également confondues :

Nouvelles représentations des données pour l'exemple numérique 2

Fig. 83 Nouvelles représentations des données pour l’exemple numérique 2

Pour ces deux exemples, l’application d’un algorithme de classification comme k-means sur les nouvelles représentations des données afin d’obtenir 2 groupes devrait donner bien plus souvent les bons résultats.

Comment obtenir un graphe de similarités

Examinons maintenant la première étape de l’algorithme de classification spectrale, c’est à dire le passage d’une matrice de similarités \(\mathbf{S}\) à un graphe de similarités \(G\).

Une étape préalable souvent nécessaire est de définir une mesure de similarité adaptée aux données de \(\mathcal{E} = \{\mathbf{x}_i, 1 \leq i \leq N\}\). En effet, les données peuvent être de nature diverse : vecteurs (cas classique de données décrites par un ensemble de variables numériques), ensembles (par ex. l’ensemble des mots d’un texte), arbres (par ex. arbres d’analyse grammaticale de phrases), etc. La compréhension des données (quelle est leur nature, quelles données sont similaires et lesquelles ne le sont pas) est très utile pour choisir une « bonne » mesure de similarité. Des mesures de similarités sont souvent définies à partir de fonctions-noyau, comme la loi normale multidimensionnelle pour des données vectorielles. De manière générale, ce qui compte le plus est le comportement de la mesure de similarité pour des données très similaires car les données très dissimilaires ne se retrouveront pas dans le même groupe.

Lorsqu’on dispose d’une métrique (mesure de distance) considérée adéquate pour les données, il est possible d’employer la loi normale pour définir les similarités : \(s_{ij} = \exp \left(- \frac{d(\mathbf{x}_i, \mathbf{x}_j)^2}{2 \sigma^2}\right)\), où \(d(\mathbf{x}_i, \mathbf{x}_j)\) est la distance entre les données vectorielles \(\mathbf{x}_i\) et \(\mathbf{x}_j\). Le choix de l’écart-type \(\sigma\) est important car, la loi normale étant à décroissance rapide, une valeur trop faible produira trop de valeurs nulles pour les \(s_{ij}\), donc nécessairement un graphe de similarités trop peu connecté. L’écart-type \(\sigma\) est en général choisi proche de la distance la plus grande entre deux données voisines dans l’arbre couvrant minimal des données de \(\mathcal{E}\).

Pour ce qui concerne le passage d’une matrice de similarités \(\mathbf{S}\) à un graphe de similarités \(G\), plusieurs solutions sont possibles :

  1. Une première solution est de considérer un graphe complet (où tous les sommets sont connectés) et d’affecter directement les valeurs de similarités issues de la matrice de similarités comme poids des arêtes du graphe, \(w_{ij} = s_{ij}\). Cela revient à compter exclusivement sur la mesure de similarité pour partitionner le graphe (les observations). Il faut également noter que plus le graphe est connecté, plus le coût d’exécution des étapes ultérieures de l’algorithme de classification spectrale est élevé.

  2. Une autre solution est de limiter la connectivité dans le graphe au \(\epsilon\)-voisinage des sommets : deux sommets \(v_i, v_j \in V\) sont connectés si leur similarité est supérieurs à un seuil, \(s_{ij} \geq \epsilon\). L’utilisation d’un seuil \(\epsilon\) est toutefois inadaptée aux situations où des groupes différents de données ont des densités très différentes : une valeur \(\epsilon\) faible limitera fortement la connectivité à l’intérieur des groupes peu denses, alors qu’une valeur plus élevée risque de connecter plusieurs groupes denses entre eux. Aussi, de nombreux sommets peuvent se retrouver isolés (voir la figure suivante) et le restent lors des étapes suivantes de l’algorithme.

  3. Limiter la connectivité dans le graphe aux \(k\) plus proches voisins (\(k\)ppv, k-nearest neighbors, kNN) des sommets répond à ce problème de présence de groupes de densités très différentes. Deux sommets \(v_i, v_j \in V\) sont ainsi connectés si \(v_i\) fait partie des \(k\)ppv de \(v_j\) ou \(v_j\) fait partie des \(k\)ppv de \(v_i\). Cette solution, conseillée en première approche, est relativement robuste au choix de la valeur de \(k\). En revanche, elle mène souvent à la création d’arêtes entre groupes de densités différentes, comme le montre la figure suivante.

  4. L’utilisation des \(k\) plus proches voisins mutuels (mutual k-nearest neighbors) permet de réduire le nombre d’arêtes entre groupes de densités différentes. Deux sommets \(v_i, v_j \in V\) sont connectés si chacun est dans les \(k\) plus proches voisins de l’autre : \(v_i\) fait partie des \(k\)ppv de \(v_j\) et \(v_j\) fait partie des \(k\)ppv de \(v_i\).

Alternatives pour la construction du graphe

Fig. 84 Alternatives pour la construction du graphe (exemple issu de [Lux07])

Pour comprendre l’idée directrice du choix des paramètres pour ces solutions alternatives il est important de se rappeler que la multiplicité de la valeur propre 0 de \(\mathbf{L}\) est égale au nombre de composantes connexes de \(G\). Or, l’algorithme de classification spectrale s’intéresse aux \(k\) plus petites valeurs propres de \(\mathbf{L}_{rw}\). En conséquence, le nombre de groupes issu de l’algorithme de classification spectrale ne peut être inférieur au nombre de composantes connexes du graphe de similarités. Mieux vaut donc éviter que le graphe \(G\) possède trop de composantes connexes car cela imposerait un nombre de groupes élevé ; une connectivité plus forte est préférable.

Supposons que le graphe \(G\) est construit à partir des \(\epsilon\)-voisinages et considérons que \(\epsilon\) est un seuil supérieur de distance et non un seuil inférieur de similarité. Pour des données vectorielles \(\mathbf{x}_i \in \mathbb{R}^d\) et la distance euclidienne classique, afin d’obtenir un graphe connecté (non décomposable en plusieurs composantes connexes) il faut choisir \(\epsilon \sim \left(\frac{\log(N)}{N}\right)^d\). Pour une métrique quelconque et des données non nécessairement vectorielles, \(\epsilon\) doit être proche de la distance la plus grande entre deux données voisines dans l’arbre couvrant minimal des données de \(\mathcal{E}\).

Si le graphe est construit à partir des \(k\) plus proches voisins, le choix de \(k\) doit respecter \(k \sim \log(N) + 1\). Lorsque la construction de \(G\) est faite grâce aux \(k\) plus proches voisins mutuels, la valeur de \(k\) devrait être plus élevée que pour les \(k\) plus proches voisins simples car pour une même valeur de \(k\) ces graphes basés sur les voisins mutuels ont moins d’arêtes que les graphes basés sur les voisins simples.

Enfin, si le graphe \(G\) est construit directement avec \(w_{ij} = s_{ij}\), le réglage doit porter plutôt sur la fonction de similarité, comme nous l’avons vu plus haut pour le passage d’une métrique à une mesure de similarité à l’aide de la loi normale.

Le nombre de groupes disjoints \(\mathcal{E}_1,\ldots,\mathcal{E}_k\) est noté par \(k\), à ne pas confondre avec le nombre de voisins employés dans la construction du graphe (notation traditionnelle). Pour la classification spectrale, le choix du nombre de groupes peut se faire à l’aide d’un critère spécifique, le premier « saut » dans la croissance des valeurs propres (en partant de la plus basse), appelé en anglais eigengap.

Considérons les valeurs propres de la matrice laplacienne \(\mathbf{L}_{rw}\) triées en ordre croissant, \(0 = \lambda_1 \leq \lambda_2 \leq \ldots\) La valeur de \(k\) (nombre de groupes) est choisie telle que dans cette liste triée on trouve un saut entre \(\lambda_k\) et \(\lambda_{k+1}\).

L’exemple suivant, issu de [Lux07], montre deux cas très simples dans lesquels les données sont unidimensionnelles. Les données sont représentées sur l’abscisse, alors que l’ordonnée indique le nombre de données de même abscisse. Dans le premier cas (à gauche) quatre groupes sont clairement identifiables et le graphique correspondant des valeurs propres montre un premier saut entre \(\lambda_4\) et \(\lambda_5\), puis un second entre \(\lambda_8\) et \(\lambda_9\) ; la valeur \(k=4\) sera ainsi retenue. Dans le second cas (à droite), les groupes sont nettement moins bien séparés, en revanche un premier saut reste bien visible entre \(\lambda_4\) et \(\lambda_5\) ; la valeur \(k=4\) sera choisie ici encore.

Illustration du saut dans les valeurs propres

Fig. 85 Illustration du saut dans les valeurs propres (eigengap, exemple issu de [Lux07])



AV07

Arthur, D. and S. Vassilvitskii. K-means++: The advantages of careful seeding. In Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete Algorithms, SODA’07, pages 1027–1035, Philadelphia, PA, USA, 2007. Society for Industrial and Applied Mathematics.

Bai13

Bair, E. Semi-supervised clustering methods, Wiley Interdiscip. Rev. Comput. Stat. 2013; 5(5): 349–361.

BBS14

Bouveyron, C., C. Brunet-Saumard. Model-based clustering of high-dimensional data: A review, Computational Statistics and Data Analysis, Volume 71, March 2014, Pages 52-78, ISSN 0167-9473, http://dx.doi.org/10.1016/j.csda.2012.12.008.

HTF01

Hastie, T., R. Tibshirani, and J. Friedman. The Elements of Statistical Learning, Springer (2001), 468–469.

JMF99(1,2)

Jain, A. K., M. N. Murty, and P. J. Flynn. Data clustering: a review. ACM Comput. Surv. 31, 3 (September 1999), 264-323. DOI=http://dx.doi.org/10.1145/331499.331504.

KR87

Kaufman, L. and Rousseeuw, P.J. Clustering by means of Medoids, Statistical Data Analysis Based on the :math:`L_1` Norm and Related Methods, ed. Y. Dodge, North-Holland, 1987, pp. 405–416.

Lux07(1,2,3,4,5,6,7)

von Luxburg, U. A tutorial on spectral clustering. CoRR, abs/0711.0189, 2007.

ST10

Shamir, O. and N. Tishby. Stability and model selection in k-means clustering. Machine Learning, 80(2):213–243, 2010.

VEB10

Vinh, N.X., J. Epps, and J. Bailey. Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance. J. Mach. Learn. Res., 11:2837–2854, Dec. 2010.

WW07

Wagner, S. and D. Wagner. Comparing Clusterings – An Overview. Technical Report 2006-04, Universität Karlsruhe (TH), 2007.