TP - PageRank

Le but de cette séance est d’implémenter l’algorithme de puissance itérée pour pouvoir ensuite calculer la centralité de vecteur propre dans le cas des graphes non orientés et également dans le des graphes orientés (Page Rank). Nous utiliserons le module networkx de Python.

Pour l'ensemble des graphes que nous allons manipuler dans ce TP, il faut que les sommets soient obligatoirement numéroté de $0$ à $n-1$.

import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import time
import math

Quelques fonctions utiles

# Return the adjacency matrix of G as list of lists
def adjacencyMatrix(G):
    return nx.to_numpy_array(G, nodelist=range(len(G))) 

# Return the matrix product of A and B
def multMatrix(A,B):
    return A.dot(B)

# Return the same graph as G whose vertices are relabeled from 0 to n-1
# and the mapping (old:new) names
def relabelGraph(G):
    mapping = dict(zip(G, range(0, len(G))))
    return nx.relabel_nodes(G, mapping), mapping 

# Return original name of the vertex from new name
def getOriginalName(mapping, v):
    return [k for k,u in mapping.items() if u == v][0]

# Draw the graph G. If a centrality measure is given in C, label the vertices. 
# Different layouts can be used: neato (default), fdp, patchwork, circo, ... 
# List of layout : https://graphviz.org/docs/layouts/
def drawGraph(G,C=None,prog="neato"):
    pos=nx.nx_agraph.pygraphviz_layout(G,prog=prog) # Position of vertices
    node_size=[300 for _ in G] # default size of nodes
    with_labels=True
    if C is not None: # if a centrality measure is given, label the vertices
        C = C.round(6)
        with_labels=False
        labels_nodes = {v:C[v][0] for v in G}
        max_C=max(C)[0]
        node_size=[math.exp(labels_nodes[v]*7/max_C) for v in labels_nodes] # size the vertices depending on the centrality measure
        nx.draw_networkx_labels(G,pos,labels_nodes,font_size=12,font_color="black") # draw labels
    nx.draw(G,pos=pos,node_size=node_size,with_labels = with_labels) # draw graph
    plt.show() 

Quelques graphes pour faire des tests

1. Un premier exemple de graphe non orienté à 4 sommets :

G1=nx.Graph()
G1.add_edge(0,1); G1.add_edge(3,0); G1.add_edge(1,2); G1.add_edge(2,0)
print("Matrice d'adjacence du graphe G1: \n", adjacencyMatrix(G1))
drawGraph(G1)

2. Un deuxième exemple de graphe généré aléatoirement sur le modèle Erdös-Rényi :

Rappel : Le graphe aléatoire $G_{n,p}$ construit sur le modèle Erdos-Rényi est un graphe à $n$ sommets dans lequel il existe une arête entre chaque paire de sommets avec une probabilité $p$.

Dans l'exemple suivant : $n=20$ et $p=\frac15$.

n=20; p=1/5
oriented=False
G2 = nx.generators.random_graphs.fast_gnp_random_graph(n,p,directed=oriented)
drawGraph(G2,prog="neato")

3. Un plus gros graphe : Gcit

Il s'agit d'un graphe de citations qui contient 27770 sommets et de 352807 arcs. Chaque sommet représente un article scientifique et les arcs représentent les citations : si l'article A cite l'article B, alors il y aura un arc de A vers B.

file_name="cit-HepTh.txt"
Gcit = nx.DiGraph()
with open(file_name, 'r') as data:
        for line in data:
            p = line.split(" ")       
            Gcit.add_edge(int(p[0]),int(p[1]))
Gcit,mappingGcit = relabelGraph(Gcit) 

Exercice 1 : Implémentation de l'algorithme de puissance itérée

1. Complétez la fonction powerIteration ci-dessous. Cette fonction implémente l’algorithme de puissance itérée ; elle prend en paramètre une matrice $A$, un nombre d'itérations maximum et un réel epsilon permettant de décider de l’arrêt. Cette fonction retourne un vecteur-colonne.
   À chaque étape de calcul (y compris à l'étape initiale), faites en sorte que la somme des valeurs du vecteur-colonne soit égale à 1.
   Pour info : un vecteur-colonne peut être obtenu avec np.ndarray(shape=(n,1)) où $n$ est la taille du vecteur.

def powerIteration(A, max_iterations=100, epsilon=1e-6):
    n = A.shape[0] # nombre de lignes de A = nombre de sommets du graphe

    # 1. Initialiser un vecteur-colonne C de taille n avec des valeurs uniformes (1/n)
    C = # TODO

    for i in range(max_iterations):
        # 2. Conserver une copie du vecteur C de l'itération précédente.
        C_prec = # TODO

        # 3. Appliquer une itération de la méthode de la puissance
        C = # TODO

        # 4. Normaliser le vecteur C pour que la somme de ses composantes soit égale à 1.
        # TODO
        
        # 5. Calculer la différence maximale en valeur absolue entre C et C_prec.
        # Si cette différence est inférieure à epsilon, arrêter la boucle.
        if # TODO :
            return C
        
    return C

2. Normalisation : Pourquoi l'étape de normalisation à chaque itération est-elle cruciale ? Que risquerait-il de se passer si on ne normalisait pas, ou seulement à la fin ?

TODO

3. Critère d'arrêt : Vous avez utilisé la plus grande différence en valeur absolue entre les composantes de C et C_prec. Une autre approche consiste à utiliser la somme des différences.
   Selon vous, lequel de ces deux critères est le plus strict ? Justifiez votre réponse.

TODO

4. Test : Testez votre fonction sur une matrice simple comme A = np.array([[0.5, 0.5], [0.5, 0.5]]). Quel est le vecteur propre attendu ? Votre fonction le trouve-t-elle ?

# TODO

Exercice 2 : Centralité de vecteur propre (cas non orienté)

1. Testez : Utilisez la fonction powerIteration écrite à l'exercice 1 pour calculer la centralité de vecteur propre dans le cas non orienté (utilisez le graphe G1 défini plus haut). Vous fixerez un réel epsilon permettant d’indiquer la précision du calcul. Faites appel à la fonction drawGraph pour voir le résultat de vos calculs s’afficher sur les sommets du graphe.

# TODO

2. Expérimentation : Modifiez votre fonction powerItération pour qu'elle retourne aussi le nombre d'itérations. Faites varier epsilon (1e-3, 1e-6, 1e-9) et observez l'impact sur le nombre d'itérations. Quel compromis cela met-il en évidence ?

# TODO

Exercice 3 : Page Rank (cas orienté)

1. Écrire une fonction matrixPageRank permettant de construire la matrice de passage du Page Rank qui nous permettra de calculer la centralité de vecteur propre dans le cas orienté . Cette fonction prend en paramètre un graphe orienté et un réel alpha correspondant au facteur de zap (dumping factor).

def matrixPageRank(G,alpha=0.85):
  # TODO: Implémenter la construction de la matrice P comme vu en cours
    pass  

2. Utilisez la fonction matrixPageRank pour calculer la matrice de passage du graphe disponible sur Moodle (que vous devez construire auparavant), puis calculez la valeur du Page Rank de chaque sommet à l'aide de la fonction powerIteration.

La valeur du Page Rank de chaque sommet est la suivante (arrondie à 2 décimales) :

[[0.33] [0.32] [0.31] [0.04]]

#TODO

3. Testez votre fonction sur le graphe de citations scientifiques GCit. Quel noeud est le plus important ?

#TODO

4. Rôle de $\alpha$ : Recalculez le PageRank de $G3$ pour $\alpha = 1$ et $\alpha = 0$. Que se passe-t-il pour $\alpha = 1$ ? Le résultat est-il celui attendu ? Que représente le PageRank quand $\alpha = 0$ ? Le résultat est-il logique ?

#TODO

Exercice 4

La construction de la matrice de passage $P$ dans matrixPageRank est coûteuse, notamment à cause des produits matriciels. L'objectif est de construire cette matrice $P$ de manière plus efficace, en travaillant ligne par ligne.

Pour un sommet $v$ donné, la ligne $v$ de la matrice $P$ dépend de deux cas :

a. Si $v$ n'est pas un puits ($d^+(v) > 0$) : la ligne $v$ de $P$ est $\alpha * (\frac{1}{d^+(v)}) A_v + \frac{(1-\alpha)}{n} u'$, où $A_v$ est la ligne $v$ de la matrice d'adjacence et $u'$ un vecteur ligne de 1.

b. Si $v$ est un puits ($d^+(v) = 0$) : la ligne $v$ de $P$ est $\frac{u'}{n}$.

1. Écrivez une nouvelle fonction matrixPageRankFast(G, alpha=0.85) qui implémente cette approche.

def matrixPageRankFast(G,alpha=0.85):
    # TODO: Implémenter la construction optimisée de la matrice P
    pass

2. Validation : Utilisez np.allclose(P_slow, P_fast) pour vérifier que votre nouvelle fonction produit bien la même matrice que matrixPageRank sur un petit graphe comme $G3$.

#TODO

3. Performance : Utilisez le module time pour comparer le temps d'exécution de matrixPageRank(Gcit) et matrixPageRankFast(Gcit). Calculez le gain de performance (ex: "la version rapide est X fois plus rapide").

#TODO