K-means challenge

Introduction

Le but de ce challenge est d’implémenter from scratch un alogrithme dit des k-means. Kesako ? C’est un algorithme de classification non supervisée… tout simplement. Hein quoi ??? Classification : nous allons tenter de regrouper nos n observations de \(\mathbb{R}^p\) au sein de K urnes en faisant en sorte que chaque urne soient plus ou moins homogènes, i.e., les observations d’une urne données possèdent les mêmes caractéristiques que celles des autres urnes n’ont pas. Non supervisé car nous ne connaissons pas a priori dans quelle urne affecter chaque observation.

OK j’ai compris mais qu’est ce donc le k-means. C’est donc un algorithme (très simple) de classification supervisée qui va regrouper les observations au sein d’une même urne selon la distance euclidienne usuelle. On peut généraliser d’ailleur en prenant \(\ell^1\) par exemple et l’on obtient les k-médianes. Bref passons… et voyons l’algorithme.

L’algorithme

Il est très simple et supporte plutôt bien la montée en dimension, i.e., \(n\) et/ou \(p\) devien(ne)t grands. Pour implémenter l’algorithme il nous faut donc un jeu de données, i.e., une matrice de réels de taille \(n \times p\), et \(2 \leq K < \infty\) le nombre d’urnes / de classes.

Faire une initialisation des centroides selon une méthodologie donnée (cf. plus loin)
Calculer la distance (euclidienne donc) entre chaque observations aux différents centroides
Affecter à chaque observation l’urne correspondant au centroide le plus proche
Mettre à jour les centroides
Revenir à l’étape 1 jusqu’à convergence.
Renvoyer la répartition des urnes.

En passant, vous l’aurez probablement deviné par vous même mais ce que l’on nomme ici un centroide n’est rien d’autre que la moyenne des observations contenues dans une urne.

L’initialisation

L’initialisation de l’algorithme kmeans est primordiale car nous pouvons converger vers un extrema local. Il faut donc faire attention. Plusieurs méthodes existes. La plus simple, appelée méthode de Forgy, consiste à définir les \(K\) centroides initiaux en piochant aléatoirement \(K\) observations parmi les \(n\) disponibles. D’autres approches sont également possible. Mais pour le moment restons en là !

Un cas d’étude et la fonction R qui fait déjà le job

Nous allons appliquer notre fonction sur un jeu de données standard : les données d’iris de Fisher. Ces données correspondent aux 150 longueurs et largeurs des pétales et sépales de trois variétés d’iris (les fleurs pas les yeux !). Elles sont disponibles sous R directement

data(iris)
summary(iris)

##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
##

Et la classification via l’algorithme kmeans sous R se fait simplement par

(classif <- kmeans(iris[,1:4], centers = 3))

## K-means clustering with 3 clusters of sizes 50, 38, 62
## 
## Cluster means:
##   Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1     5.006000    3.428000     1.462000    0.246000
## 2     6.850000    3.073684     5.742105    2.071053
## 3     5.901613    2.748387     4.393548    1.433871
## 
## Clustering vector:
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
##  [71] 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 2 2 2
## [106] 2 3 2 2 2 2 2 2 3 3 2 2 2 2 3 2 3 2 3 2 2 3 3 2 2 2 2 2 3 2 2 2 2 3 2
## [141] 2 2 3 2 2 2 3 2 2 3
## 
## Within cluster sum of squares by cluster:
## [1] 15.15100 23.87947 39.82097
##  (between_SS / total_SS =  88.4 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"

Visualisons le résultat de cette classification et comparons là à la vérité. Le suspens est insoutenable ;-)

par(mfcol = c(2, 6), mar = c(4, 5, 0.5, 0.5), ps = 18)
for (i in 1:3){
  for (j in (i+1):4){
    plot(iris[,i], iris[,j], col = classif$cluster, xlab = names(iris)[i],
         ylab = names(iris)[j])
    points(classif$centers, pch = 15, col = 1:3)
    
    plot(iris[,i], iris[,j], col = as.numeric(iris[,5]), xlab = names(iris)[i],
         ylab = names(iris)[j])
  }
}

Plutôt pas mal non ?

Application sur MNIST

Le jeu de données MNIST est très connu dans le domaine de l’apprentissage statistique. Ici nous allons voir comment la classification par kmeans se débrouille.