6  Dbscan et clustering spectral

6.1 L’algorithme DBSCAN

L’algorithme DBSCAN (Density Based Spatial Clustering of Applications with Noise, Ester et al. (1996)) fait partie des méthodes basées sur la densité : les clusters correspondent à des zones de fortes densité séparées par des zones où la densité est plus faible. Ces zones sont définies par deux types de points :

  • les noyaux : des points qui contiennent plus de minPts points à une distance inférieure à epsilon ;
  • les points de bordure : les points qui vont se situés en bordure des clusters ou qui sont isolés.

On charge le package

library(dbscan)

Exercice 6.1 (Noyaux et points de bordure) On considère le “nuage” de points suivant :

set.seed(1234)
tbl <- tibble(X1=c(rnorm(5,mean=0.25,sd=0.05),0.75,0.5),
              X2=c(rnorm(5,mean=0.45,sd=0.05),0.4,0.2),
              label=as.character(1:7))

              
ggplot(tbl)+aes(x=X1,y=X2)+geom_point()+
  geom_text(aes(label=label),vjust=2)

  1. On fixe eps=0.13 et minPts=4. À l’aide de calculs simples, identifier les noyaux et points de bordure.

    Il suffit de calculer les distances entre individus et de compter le nombre de points à une distance inférieurs de eps de chaque individu.

    dist(tbl[,1:2])
               1          2          3          4          5          6
2 0.09181292                                                       
3 0.12608683 0.04037506                                            
4 0.07814202 0.13115736 0.17150926                                 
5 0.10754298 0.01749459 0.03699214 0.13969293                      
6 0.56539041 0.48659333 0.44635393 0.61766891 0.47857535           
7 0.41486238 0.32359543 0.29649658 0.42904942 0.30734626 0.32015621

    On déduit que 1, 2, 3 et 5 sont des noyaux, 4, 6 et 7 sont des points de bordure.

  2. Retrouver ces résultats à l’aide de la fonction is.corepoint.

    is.corepoint(tbl[,1:2],eps = 0.13,minPts = 4)
    [1]  TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE

  3. Effectuer l’algorithme dbscan avec ces valeurs de paramètre et interpréter.

    (db <- dbscan(tbl[,1:2],eps = 0.13,minPts = 4))
    DBSCAN clustering for 7 objects.
Parameters: eps = 0.13, minPts = 4
Using euclidean distances and borderpoints = TRUE
The clustering contains 1 cluster(s) and 2 noise points.

0 1 
2 5 

Available fields: cluster, eps, minPts, dist, borderPoints
    db$cluster
    [1] 1 1 1 1 1 0 0

    L’algorithme a identifié 1 cluster de 5 points et 2 outliers (6 et 7).

  4. Est-ce que des points de bordure peuvent être affectés à des cluster ? Justifier.

    Oui ! On voit par exemple que 4 est dans le cluster. Ce point est en effet connecté aux 4 autres. On voit par exemple qu’il est connecté avec 2 car 4 et 2 sont tous les deux accessibles depuis 1.

Exercice 6.2 (Calibration de dbscan) On reprend les données de l’Exercice 5.1 :

tbl <- read_delim("data/donclassif.txt",delim = ";")
ggplot(tbl)+aes(x=V1,y=V2)+geom_point()

En utilisant la stratégie proposée dans l’aide de dbscan, calibrer l’algorithme pour essayer d’identifier au mieux les différents clusters. On pourra envisager une deuxième étape pour affecter les petits clusters aux gros…

On fixe minPts=3 et on regarde la distance entre chaque point et son deuxième plus proche voisin.

kNNdistplot(tbl,k=2)

On choisit eps=0.3 et on visualise les résultats :

(db1 <- dbscan(tbl,eps=0.3,minPts = 3))
DBSCAN clustering for 2800 objects.
Parameters: eps = 0.3, minPts = 3
Using euclidean distances and borderpoints = TRUE
The clustering contains 12 cluster(s) and 39 noise points.

  0   1   2   3   4   5   6   7   8   9  10  11  12 
 39 400 400 400 400 400 398 338   4   8   3   7   3 

Available fields: cluster, eps, minPts, dist, borderPoints
noyau <- is.corepoint(tbl,eps=0.3,minPts = 3)
tbl_db <- tbl |> mutate(dbscan=as.factor(db1$cluster),noyau=noyau)
ggplot(tbl_db)+aes(x=V1,y=V2,color=dbscan,shape=noyau)+geom_point()

On remarque que les groupes 8 à 12 ont des effectifs très faibles. On propose de les affecter au “gros” groupe le plus proche en utilisant la distance du min :

library(LearnClust)
D <- rep(0,7)
tbl_db <- tbl_db |> mutate(new=dbscan)
for (i in c(0,8:12)){
  tbl1 <- tbl_db |> filter(dbscan %in% i) |> select(1:2) |> as.matrix()
  for (j in 1:7){
    tbl2 <- tbl_db |> filter(dbscan %in% j) |> select(1:2) |> as.matrix()
    D[j] <- clusterDistance(tbl1,tbl2,'MIN','MAN')
    tbl_db$new[tbl_db$dbscan==i] <- which.min(D)
  }
}
ggplot(tbl_db)+aes(x=V1,y=V2,color=new)+geom_point()

6.2 Clustering spectral

Le clustering spectral est un algorithme de classification non supervisé qui permet de définir des clusters de nœuds sur des graphes ou d’individus pour des données individus/variables. L’algorithme est basé sur la décomposition spectrale du Laplacien (normalisé) d’une matrice de similarité, il est résumé ci-dessous :

Entrées :

  • tableau de données \(n\times p\)
  • \(K\) un noyau
  • \(k\) le nombre de clusters.
  1. Calculer la matrice de similarités \(W\) sur les données en utilisant le noyau \(K\)
  2. Calculer le Laplacien normalisé \(L_{\text{norm}}\) à partir de \(W\).
  3. Calculer les \(k\) premiers vecteurs propres \(u_1,\dots,u_k\) de \(L_{\text{norm}}\). On note \(U\) la matrice \(n\times k\) qui les contient.
  4. Calculer la matrice \(T\) en normalisant les lignes de \(U\) : \(t_{ij}=u_{ij}/(\sum_\ell u_{i\ell}^2)^{1/2}\).
  5. Faire un \(k\)-means avec les points \(y_i,i=1,\dots,n\) (i-\(ème\) ligne de \(T\)) \(\Longrightarrow\) \(A_1,\dots,A_k\).

Sortie : clusters \(C_1,\dots,C_k\) avec \[C_j=\{i|y_i\in A_j\}.\]

L’objet de ce chapitre est de travailler sur cet algorithme en le programmant, puis en utilisant la fonction specc du package kernlab.

On crée tout d’abord un graphe avec trois composantes connexes : on utilise la commande sample_gnp() qui permet de créer un graphe selon le modèle d’Erdos-Renyi.

library(igraph)
set.seed(1)
n1 <- 5
n2 <- 3
n3 <- 2
n <- n1+n2+n3
# il faut prendre des grandes valeurs de p sinon on risque d'avoir des sous-graphes non connexes 
p1 <- 0.85
p2 <- 0.75
p3 <- 0.7
G1 <- sample_gnp(n1,p1)
G2 <- sample_gnp(n2,p2)
G3 <- sample_gnp(n3,p3)
G <- G1 + G2 + G3 # il cree un graphe avec ces 3 sous-graphes
plot(G)

On vérifie le nombre de composantes connexes

components(G)$no
[1] 3

Exercice 6.3 (Laplacien non normalisé)

  1. Calculer la matrice d’adjacence de G et en déduire le Laplacien normalisé. On pourra utiliser la fonction as_adj.

    A <- as_adj(G,sparse=F)
D <- diag(rowSums(A))
D_moins1_2 <- diag(1/sqrt(diag(D))) 
LN <- diag(n) - D_moins1_2 %*% A %*%D_moins1_2

  2. Retrouver ce Laplacien avec la fonction laplacian_matrix.

    LN <- laplacian_matrix(G,norm=TRUE,sparse=F)

  3. Calculer les valeurs propres et représenter les sur un graphe. Que remarquez-vous ?

    specN <- eigen(LN)
specN$values
     [1] 2.000000e+00 2.000000e+00 1.564333e+00 1.333333e+00 1.250000e+00
 [6] 1.000000e+00 8.523332e-01 1.110223e-15 1.110223e-15 8.881784e-16

    On observe un trou spectral entre les valeurs propres 3 et 4. Conformément à la théorie, l’ordre de multiplicité de la valeur propre 0 est égal au nombre de composantes connexes du graphe.

    dfN <- tibble(vp=1:length(specN$values),valeur=rev(specN$values))
ggplot(dfN)+aes(x=vp,y=valeur)+geom_point()+theme_classic()

  4. Obtenir les trois vecteurs propres associés à la valeur propre nulle. Commenter.

    On calcule les 3 vecteurs propres :

    U <- specN$vectors[,n:(n-2)]
U
                [,1]      [,2]      [,3]
 [1,] -0.5000000 0.0000000 0.0000000
 [2,] -0.4330127 0.0000000 0.0000000
 [3,] -0.3535534 0.0000000 0.0000000
 [4,] -0.5000000 0.0000000 0.0000000
 [5,] -0.4330127 0.0000000 0.0000000
 [6,]  0.0000000 0.7071068 0.0000000
 [7,]  0.0000000 0.5000000 0.0000000
 [8,]  0.0000000 0.5000000 0.0000000
 [9,]  0.0000000 0.0000000 0.7071068
[10,]  0.0000000 0.0000000 0.7071068

  5. Normaliser ces vecteurs. On pourra utiliser la fonction

    normalize <- function(x){
  return(x/sqrt(sum(x^2))) 
}

    On normalise à l’aide de la fonction normalize :

    U.norm <- t(apply(U,1,normalize)) 
U.norm
          [,1] [,2] [,3]
 [1,]   -1    0    0
 [2,]   -1    0    0
 [3,]   -1    0    0
 [4,]   -1    0    0
 [5,]   -1    0    0
 [6,]    0    1    0
 [7,]    0    1    0
 [8,]    0    1    0
 [9,]    0    0    1
[10,]    0    0    1

  6. Terminer l’algorithme avec le \(k\)-means.

    res <- kmeans(U.norm,3,nstart=100)
res$cluster
     [1] 1 1 1 1 1 2 2 2 3 3

Exercice 6.4 (Clustering spectral avec specc) On reprend les données de l’Exercice 5.1 :

tbl <- read_delim("data/donclassif.txt",delim = ";")

Faire le “meilleur” clustering spectral, on pourra utiliser plusieurs noyaux et plusieurs valeurs de paramètres.

  1. On commencera par lancer la procédure sur un sous-échantillon de taille 1000

  2. On proposera ensuite une autre solution qui permettra de traiter tous les observations.

On commence en tirant un sous-échantillon et un choisit 9 groupes.

set.seed(1234)
perm <- sample(nrow(tbl))
tbl1 <- tbl[perm[1:800],]
library(kernlab)
g1 <- specc(as.matrix(tbl1),centers=9,kpar=list(sigma=100))
tbl2 <- tbl1 |> mutate(g1=as.factor(g1))
ggplot(tbl2)+aes(x=V1,y=V2,color=g1)+geom_point()

Pour pouvoir appliquer l’algorithme à toutes les données on

  1. fait un \(k\)-means avec 1000 groupes

  2. fait le clustering spectral sur les centroïdes

  3. affecte tous les individus en fonction du \(k\)-means initial.

e1 <- kmeans(tbl,centers=1000,nstart = 100,iter.max = 100)
centres <- e1$centers
g1 <- e1$cluster
e2 <- specc(centres,centers=9,kernel="rbfdot",
            kpar=list(sigma=100))
g2 <- e2[g1]
tbl1 <- tbl |> mutate(KM=g1,final=as.factor(g2))
ggplot(tbl1)+aes(x=V1,y=V2,color=final)+geom_point()