Le projet consiste à implémenter en R l’algorithme MCL de partitionnement de graphe (clustering) : à partir d’un graphe et d’un paramètre nommé inflate factor (noté IF par la suite), l’algorithme découpe le graphe en clusters. Exemple sur le graphe ci-dessous avec un IF de 3 :

Introduction sur les graphes

Un graphe est un objet mathématique composé de sommets (vertex) reliés par des arètes (edge). Une des représentations informatiques consiste en une matrice carrée nommée matrice d’adjacence (adjacency matrix). Une cellule de cette matrice contient le poids de l’arète reliant deux sommets ou zéro en absence de lien.

Exemple :

    A  B  C  D
A   0  1  0  4
B   1  0  1  1
C   0  1  0  0
D   1  4  0  0

correspond au graphe :

Remarque : Dans ce dessin, l’épaisseur du trait de l’arète est proportionnel à son poids.

On distingue les graphes orientés (directed graph) des graphes non orientés (undirected graph). Dans un graphe orienté, les arètes sont appelées arcs et ont une direction (de A vers B). La matrice d’adjacence n’est donc plus obligatoirement symétrique et un choix arbitraire est fait concernant son interprétation : par exemple, les sommets sources correspondent aux lignes et les sommets destinations correspondent aux colonnes. Illustration :

    A  B  C  D                                   A  B  C  D
A   0  1  0  4                               A      1     4
B   0  0  0  2   ou de manière plus lisible  B            2
C   0  1  1  0                               C      1  1  
D   0  0  0  0                               D

Remarque : Le lien de C vers lui-même s’appelle une boucle (loop).

Principes et algorithme de MCl

Une bonne méthode de clustering maximise

  • la similarité des objets appartenant à un même cluster (cohésion)
  • la dissemblance entre objets appartenant à des clusters différents (séparation).

En ce qui concerne le partitionnement de graphes, une bonne méthode vise à ce que les clusters obtenus

  • aient un maximum d’arètes ou d’arcs reliant les sommets d’un même cluster
  • aient un minimum d’arètes ou d’arcs reliant les sommets de clusters différents

Pour identifier les clusters, MCL se base sur la remarque suivante : si on place un marcheur aléatoire sur un sommet d’un cluster, il a plus de chances (en se promenant aléatoirement d’un sommet à un autre) de rester dans le cluster (plus grande densité de liens) que d’en sortir.

Le procédé s’appuie donc sur le calcul des probabilités de marches aléatoires (probabilité de passer par un sommet et probabilité d’emprunter un arc) dans un graphe donné. Le calcul s’effectue sur des matrices de Markov qui représentent les probabilités de transition d’un sommet à un autre ; il s’agit donc d’une matrice d’adjacence dont la somme des lignes vaut 1 (la somme des probabilités de sortie d’un sommet vaut 1). Exemple avec le graphe précédent :

Matrice d'adjacence                Matrice de Markov
    A  B  C  D                        A   B    C    D
A   0  1  0  4                    A   0  0.2   0   0.8
B   0  0  0  2                    B   0   0    0    1
C   0  1  1  0                    C   0  0.5  0.5   0
D   0  0  0  0                    D   0   0    0    0

Remarque : Ces matrices sont également appelées matrices stochastiques ou encore matrices de transition.

L’algorithme MCL simule des marches aléatoires dans un graphe en alternant deux oprations appelées expansion et inflation. L’expansion correspond au calcul de marches aléatoires de grandes longueurs et coincide à élever une matrice stochastique à une certaine puissance. Concrètement, l’expansion consiste à multiplier la matrice avec elle-même avec le produit matriciel. L’inflation consiste à exagérer les probabilités de marches au sein d’un même cluster et à atténuer les probabilités de marches inter-clusters. En pratique, l’inflation consiste à élever chaque cellule de la matrice à une certaine puissance (inflate factor) puis à normaliser les valeurs obtenues afin que la matrice soit stochastique.

Au final, à force de répéter les opérations d’expansion et d’inflation sur la matrice de transition (et donc sur le graphe), celui-ci est décomposé en différentes composantes (composantes connexes) déconnectées les unes des autres et qui correspondent aux clusters. En d’autres termes, l’algorithme converge et la matrice n’évolue plus.

L’algorithme est donc le suivant :

M : matrice d’adjacence

ajouter les boucles à M
M_1 la matrice stochastique obtenue à partir de M
tant que on observe un changement dans la matrice
faire :
   M_2 \(\leftarrow\) expansion(M_1)
   M_2 \(\leftarrow\) inflation(M_2, inflate_factor)
   déterminer s’il y a eu un changement
   M_1 \(\leftarrow\) M_2

Le résultat est la matrice d’adjacence M_1 ou bien les composantes connexes du graphe correspondant à cette matrice.

Mise en oeuvre en R

Prise en main de la librairie igraph

Pour le chargement, le dessin et l’analyse des graphes, vous pourrez utiliser la librairie igraph.

Pour installer la librairie (si nécessaire) :

install.packages('igraph')

Puis pour l’utiliser :

library(igraph)

Pour créer un graphe à partir d’une matrice d’adjacence :

Exemple pour un graphe non orienté :

Matrice d'adjacence:
    A  B  C  D
A   0  1  0  4
B   1  0  1  1
C   0  1  0  0
D   1  4  0  0

Création de la matrice :

M = matrix(c(0,1,0,1, 1,0,1,4, 0,1,0,0, 4,1,0,0), ncol=4)
colnames(M) = rownames(M) = c('A', 'B', 'C', 'D')
M
##   A B C D
## A 0 1 0 4
## B 1 0 1 1
## C 0 1 0 0
## D 1 4 0 0

Création du graphe à partir de la matrice

g=graph.adjacency(M, mode='undirected', weighted=TRUE)

La fonction V(g) permet d’accéder aux sommets (vertices) et la fonction E(g) aux arêtes (edges). On peut ainsi attribuer une étiquette aux sommets :

V(g)$label=c('A','B','C','D')
V(g)
## + 4/4 vertices, named, from 33cee8f:
## [1] A B C D
V(g)$label
## [1] "A" "B" "C" "D"

On peut accéder au poids des arêtes :

E(g)
## + 4/4 edges from 33cee8f (vertex names):
## [1] A--B A--D B--C B--D
E(g)$weight
## [1] 1 4 1 4

Remarque : il est possible de la même manière d’ajouter n’importe quel attribut sur les sommets ou les arêtes du graphe.

Et tracer le graphe :

plot(g, vertex.color='white', edge.width=E(g)$weight)

Pour un graphe orienté :

    A  B  C  D
A   0  1  0  4 
B   0  0  0  2  
C   0  1  1  0
D   0  0  0  0

Création de la matrice :

M = matrix(c(0,0,0,0, 1,0,1,0, 0,0,1,0, 4,2,0,0), ncol=4)

Création du graphe et ajout des étiquettes :

g=graph.adjacency(M, mode='directed', weighted=TRUE)
V(g)$label=c('A','B','C','D')

Dessin du graphe :

plot(g, vertex.color='white', edge.width=E(g)$weight)

Chargement du graphe fourni pour le projet : celui-ci est dans un format très simple appelé edgelist dans igraph. Sur chaque ligne du fichier, une arête est décrite par le sommet source et le sommet destination :

g = read.graph('http://silico.biotoul.fr/enseignement/m1/math/projet/toy.tgr', format='ncol', directed=FALSE)

Pour afficher le graphe :

plot(g)

Pour obtenir la matrice d’adjacence du graphe :

m = get.adjacency(g)
m
## 37 x 37 sparse Matrix of class "dgCMatrix"
##    [[ suppressing 37 column names '0', '1', '2' ... ]]
##                                                                             
## 0  . 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 . . . . . . . . . . . . . . . . .
## 2  1 1 . . 1 . . . . . . . . 1 . . 1 . 1 . 1 1 1 1 1 . . . . . . . . . . . .
## 3  1 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 4  1 1 1 . . . 1 1 1 1 1 1 1 1 1 1 1 1 . 1 . . . . . 1 . . . . . . . . . . .
## 5  1 1 . . . . 1 . . 1 1 . 1 1 1 1 . 1 . . . . . . . . . . . . . . . . . . .
## 6  1 1 . . 1 1 . . . 1 1 . 1 1 1 1 . 1 . . . . . . . . 1 . . . . . . . . . .
## 7  1 1 . . 1 . . . . . 1 . . 1 1 . 1 . . 1 . . . . . . . 1 1 1 1 1 1 . . . .
## 8  1 . . . 1 . . . . . 1 1 . . 1 . . . . . . . . . . . . . . . . . . . . . .
## 9  1 1 . . 1 1 1 . . . 1 . 1 1 1 1 . 1 . . . . . . . . 1 . . . . . . . . . .
## 10 1 1 . . 1 1 1 1 1 1 . 1 1 1 1 1 1 1 . 1 . . . . . 1 . . . . . . . 1 . . .
## 11 1 1 . . 1 . . . 1 . 1 . . 1 . . . . . . . . . . . . . . . . . . . . . . .
## 12 1 1 . . 1 1 1 . . 1 1 . . 1 1 1 . 1 . . . . . . . . 1 . . . . . . . . . .
## 13 1 1 1 . 1 1 1 1 . 1 1 1 1 . 1 1 1 1 . . . . . . . . . . . . . . . . . . .
## 14 1 1 . . 1 1 1 1 1 1 1 . 1 1 . 1 . 1 . 1 . . . . . . . . . . . . . 1 . . .
## 15 1 1 . . 1 1 1 . . 1 1 . 1 1 1 . . 1 . . . . . . . . . . . . . . . . . . .
## 16 1 1 1 . 1 . . 1 . . 1 . . 1 . . . . 1 . 1 1 . 1 1 . . . . . . . . . . . .
## 17 1 1 . . 1 1 1 . . 1 1 . 1 1 1 1 . . 1 . . . . . 1 . . . . . . . . . . . .
## 18 . 1 1 . . . . . . . . . . . . . 1 1 . . 1 1 1 1 1 . . . . . . . . . . . .
## 19 . 1 . . 1 . . 1 . . 1 . . . 1 . . . . . . . . . . . . 1 1 1 1 1 1 . . . .
## 20 . . 1 . . . . . . . . . . . . . 1 . 1 . . 1 1 1 1 . . . . . . . . . . . .
## 21 . . 1 . . . . . . . . . . . . . 1 . 1 . 1 . . . . . . . . . . . . . . . .
## 22 . . 1 . . . . . . . . . . . . . . . 1 . 1 . . 1 1 . . . . . . . . . . . 1
## 23 . . 1 . . . . . . . . . . . . . 1 . 1 . 1 . 1 . 1 . . . . . . . . . . . 1
## 24 . . 1 . . . . . . . . . . . . . 1 1 1 . 1 . 1 1 . . . . . . . . . . . . 1
## 25 . . . . 1 . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . .
## 26 . . . . . . 1 . . 1 . . 1 . . . . . . . . . . . . . . . . . . . . . . . .
## 27 . . . . . . . 1 . . . . . . . . . . . 1 . . . . . . . . 1 1 1 1 1 . . . .
## 28 . . . . . . . 1 . . . . . . . . . . . 1 . . . . . . . 1 . 1 1 1 1 . . . .
## 29 . . . . . . . 1 . . . . . . . . . . . 1 . . . . . . . 1 1 . 1 1 1 . . . .
## 30 . . . . . . . 1 . . . . . . . . . . . 1 . . . . . . . 1 1 1 . 1 1 . . . .
## 31 . . . . . . . 1 . . . . . . . . . . . 1 . . . . . . . 1 1 1 1 . 1 . . . .
## 32 . . . . . . . 1 . . . . . . . . . . . 1 . . . . . . . 1 1 1 1 1 . . . . .
## 33 . . . . . . . . . . 1 . . . 1 . . . . . . . . . . . . . . . . . . . 1 1 .
## 34 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . 1 .
## 35 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 . .
## 36 . . . . . . . . . . . . . . . . . . . . . . 1 1 1 . . . . . . . . . . . .

Implémentation de MCl

Vous allez donc implémenter MCL en R comme décrit plus haut. Pour cela, vous définirez plusieurs fonctions :

  • add_loops(M): qui prend en paramètre une matrice d’adjacence et qui renvoie la même matrice à laquelle sont ajoutées les boucles. La fonction supposera que la matrice d’adjacence est composée de 0 et de 1.
  • make_stochastic(M): qui prend en paramètre une matrice d’adjacence et qui renvoie la matrice stochastique correspondante
  • expansion(M): qui prend en paramètre une matrice stochastique et qui renvoie la matrice après l’opération d’expansion
  • inflation(M, inflate=2): qui prend en paramètre la matrice renvoyée par la fonction précédente et renvoie la matrice après l’opération d’inflation. Elle prend également le paramètre inflate correspondant à l’inflate factor (et qui par défaut sera de 2)
  • MCL(M, inflate=2) : qui prend en paramètre une matrice d’adjacence (avec des 0 et des 1) et l’inflate factor (par défaut 2) et qui renvoie la matrice d’adjacence (arrondie à 3 chiffres après la virgule) après application de l’algorithme de MCL

Pour la dernière fonction, vous aurez aussi besoin d’une fonction qui détecte s’il y a eu un changement. Pour cela, vous allez utilisez une fonction chaos qui calcule la valeur suivante :

  • soit max la valeur maximale sur une ligne
  • soit sum_sq la somme des valeurs de la ligne élevées au carré (chaque élément est élevé au carré puis on en fait la somme)
  • on a donc row_chaos = max - sum_sq : différence entre le maximum de la ligne et la somme des carrées des éléments pour chaque ligne.

La fonction renverra la plus forte valeur observée sur les lignes (le maximum des row_chaos obtenus sur chacune des lignes).

Remarque : Cette fonction chaos tient sur une ligne de code R.

Le pseudo-code la fonction MCL devient :

MCL(M, inflate=2) :
   M2 = add_loops(M)
   M2 = make_stochastic(M2)
   change = 1
   tant que change > 0.001 faire :
      M2 = expansion(M2)
      M2 = inflation(M2, inflate)
      change = chaos(M2)
   renvoyer arrondi(M2,3)

Ainsi, en appelant la fonction MCL comme suit :

g = read.graph('http://silico.biotoul.fr/enseignement/m1/math/projet/toy.tgr', format='ncol', directed=FALSE)

m = MCL(get.adjacency(g), inflate=3)

Création du graphe à partir de la matrice renvoyée par MCL

g2=graph.adjacency(m, mode='undirected', weighted=TRUE)

Le graphe obtenu est donc le suivant :

plot(g2)

La fonction clusters nous renvoie les composantes connexes du graphes (les clusters) :

clusters(g2)
## $membership
##  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 
##  1  1  2  1  1  1  1  3  1  1  1  1  1  1  1  1  2  1  2  3  2  2  2  2  2  1  1  3  3  3  3  3  3  4  4  4  2 
## 
## $csize
## [1] 17  9  8  3
## 
## $no
## [1] 4

Ainsi, clusters(g2)$membership est un vecteur qui à chaque sommet associe un numéro de cluster. Nous l’utilisons pour colorer les sommets du graphe de départ :

plot(g, vertex.color=clusters(g2)$membership)

Dossier à rendre

Un compte rendu scientifique (habituel) sous la forme de fichier R-Markdown avec les fichiers .Rmd et .html ; il doit contenir notamment :

  • Les fonctions réalisées et commentées
  • Le chargement et le clustering du graphe fourni (toy.tgr) avec différentes valeurs d’inflate factor avec les dessins associés
  • Le chargement et le clustering du graphe de votre projet avec différentes valeurs d’inflate factor avec les dessins associés

Questions à aborder :

  • Que se passe-t-il quand on augmente l’inflate factor ?
  • Quelle est la valeur minimale pour ce paramètre ?
  • Quelle est la valeur maximale pour ce paramètre ?
  • Quel est, selon vous, le meilleur IF pour ces graphes ? et pourquoi ?
  • Comment, selon vous, le déterminer ?

Le compte rendu doit aussi comprendre une partie bilan avec par exemple :

  • les problèmes rencontrés,
  • des perspectives sur l’utilisation d’une telle méthode

Le tout est à envoyer par mail avant les fêtes de fin d’années.

Liens