PAGERANK#

Description#

L’algorithme PageRank de Google qui lui permet de classer les pages Web par ordre de pertinence est basé sur les chaînes de Markov. Il part du principe que si un utilisateur se promène de liens en liens sur des pages comportant le mot-clé recherché, la page par laquelle il passe le plus souvent est sans doute la plus pertinente. Formellement, si $N$ pages numérotées de $1$ à $N$ réfèrent toutes à un mot-clé donné, on définit

\[\begin{split}\begin{align*} N_{ij}&= \begin{cases} 1 & \text{ s'il y a un lien de $i$ vers $j$},\\ 0 & \text{ sinon.} \end{cases}\\ N_i&=\sum_{j\neq i} N_{ij} \end{align*}\end{split}\]

On construit une chaîne de Markov, construit sur $\{1,\cdots, N\}$, de matrice de transition

\[p_{ij}=\frac{N_{ij}}{N_i}\cdotp\]

Cela revient à imaginer une marche aléatoire sur les pages concernées où l’on choisit au hasard une page parmi celles qui sont possibles sur la page courante. Comme il peut arriver que des pages soient isolées, c’est-à-dire sans liens avec les autres où que des sous-ensembles de pages soient sans contact, autrement dit que la matrice ne soit pas irréductible; on ne travaille pas directement avec la marche aléatoire sur le graphe mais sur une marche un peu “modifiée” . Soit $\Delta$ la matrice qui comporte des $0$ sur la diagonale et des $1/(N-1)$ partout ailleurs. La matrice de transition considérée par PageRank est la matrice

\[\alpha P + (1-\alpha) \Delta,\]

où $\alpha$ a été fixée empiriquement à $0{,}85$.

Votre travail#

La maison vous offre le code de génération du graphe à étudier. Vous devez impérativement utiliser cette fonction avec un nombre de noeuds égal à 15 et 150 et une graîne de génération aléatoire égale à ce que vous voulez.

Pour chacun de ces deux graphes censé représenter les connexions entre les pages associées à un mot-clé, vous devez

Calculer algébriquement la probabilité stationnaire de la marche aléatoire. Soit $\pi_m$ le vecteur correspondant.
Classer les sommets selon leur probabilité stationnaire : le premier sommet est celui qui correspond à la composante de la probabilité stationnaire la plus grande, le deuxième est associé à la deuxième plus grande composante, etc. Soit $\pi_m^\prec$ ce vecteur.
Simuler des marches aléatoires sur votre graphe de $10^n$ pour $n=3,4,5,6$, pas et calculer la fréquence des passages dans chaque sommet tout le long de votre marche. Soit $\pi_s^n$ le vecteur correspondant.
Calculer la norme euclidienne de $\pi_s^n-\pi_m$ pour $n=3,4,5,6$.
Pour chaque valeur de $n$, classer les sommets par ordre de fréquence décroissant. Soit $\pi_s^{n,\prec}$ ce vecteur.
Pour vérifier que cet ordre correspond bien au précédent malgré l’imprécision due à la simulation, on calcule le nombre de composantes où l’écart relatif entre $\pi_m^\prec$ et $\pi_s^{n,\prec}$ est supérieur à 1%.
Pour comparer les durées entre les simulation et la résolution algébrique, on prendra comme unité de temps la durée de la résolution du système.
Au final, quels sont les avantages et inconvénients de chacune des deux méthodes en termes de précision, de rapidité, de robustesse si le nombre de noeuds augmente ou de tout autre critère qui vous viendrait à l’esprit ?

Pour commencer#

La page qui s’affiche actuellement est un notebook constitué de “cellules” de deux types : des cellules de texte (comme celle-ci, en langage Markdown, utilisé par exemple par Wikipedia) et des cellules de code (en langage python). Ce notebook peut être exportée sous forme de page Web, lisible par n’importe quel navigateur, ou de document pdf (voir le menu File).

Vous pouvez naviguer entre les cellules avec les flèches Haut et Bas (essayez !), éditer chaque cellule en double-cliquant dessus, et exécuter la cellule en cours en tapant Ctrl-Enter. L’ensemble des raccourcis clavier s’affiche en cliquant sur la touche h.

Voici pour démarrer deux cellules à compléter, une de texte puis une de code :

Choisissez votre graîne#

# Quel est votre graine ?
# Celle-ci servira à générer un graphe qui vous est propre
graine=1

Chaque cellule de code fonctionne comme une console python classique. Vous pouvez en particulier utiliser la complétion (avec Tab) et demander de l’aide en ligne. Par exemple :

import numpy as np
?np.dot
?np.random

Génération du graphe#

On utilise un modèle de génération de graphes dit de Erdös et Renyi. On part d’un graphe dirigé complet à $nbVertices$ sommets et chacune des $nbVertices(nbVertices-1)$ arêtes dirigées est conservée avec probabilité $p$ prise ici à $2/nbVertices$ pour avoir un graphe dense mais pas trop.

import networkx as nx
import numpy as np


def createGraph(nbVertices,randomSeed):
    G=nx.erdos_renyi_graph(nbVertices,2./nbVertices,seed=randomSeed,directed=True)
    M=nx.adjacency_matrix(G).toarray().astype(float)
    for i in np.arange(nbVertices):
        u=np.sum(M[i,:])
        if (u>0):
            M[i,:]/=u
    U=(np.ones_like(M)-np.identity(nbVertices))/(nbVertices-1)
    Unorm=np.sum(U[0,:])
    U=U/Unorm
    M=0.85*M+0.15*U
    for i in np.arange(nbVertices):
        u=np.sum(M[i,:])
        if (u>0):
            M[i,:]/=u
    return M

Visualisation#

Pour avoir une idée de ce qu’est votre graphe, vous pouvez le visualisez grâce à la commande draw_networkx

%matplotlib inline
nbVertices=15
import matplotlib.pyplot as plt
G=nx.erdos_renyi_graph(nbVertices,2./nbVertices,seed=graine,directed=True)
plt.axis('off')
nx.draw_networkx(G,with_labels=True)

Calculs#

Créer une fonction qui calcule la probabilité stationnaire par résolution du système d’équation $\pi=\pi P$ et $\sum_{i=1}^{nbVertices}\pi(i)=1$. Cette fonction retourne le vecteur de probabilité stationnaire mais aussi l’ordre de ces composantes. Pour un vecteur $[0.2,0.4,0.3]$, on veut le vecteur $[1,2,0]$. Rappel : en Python, les tableaux sont numérotés à partir de $0$. Rien de difficile ici, regarder la documentation des fonction numpy.argsort et numpy.flipud
Créer une fonction qui simule la marche aléatoire avec un nombre de pas passé en paramètre. Cette fonction retourne le vecteur des fréquences empiriques ainsi l’ordre de ses composantes.

?np.argsort
?np.flipud

def with_eigenvals(b):
    ?np.hstack
    ?np.linalg.solve

def random_walk(b,K):
    ?np.random.choice

import time
# Crée le ficher numero_de_casier.txt
with open(str(casier)+'.txt','w') as f:
    for noeuds in [15,150]:
        M=createGraph(noeuds,casier)
        u0=time.time()
        a,b=with_eigenvals(M)
        t=time.time()-u0
        # N° casier, nb de noeuds, temps de résolution du système linéaire
        f.write("\n" +str(casier)+"\t"+str(noeuds)+"\t"+"{0:.2E}".format(t)+"\t")
        for K in [1000,10000,100000,1000000]:
            u0=time.time()
            c,d=random_walk(M,K)
            # Unité de temps = temps de résolution du système
            # Durée de la simulation
            u=(time.time()-u0)/t
            f.write("{0:.1f}".format(u)+"\t")
            # Norme de la différence entre pi_s et pi_m
            f.write("{0:.4f}".format(np.linalg.norm(b-d))+"\t")
            # Nombre de composantes ordonnées différentes (avec une erreur relative supérieure à 1%)
            f.write(str(np.sum(np.abs((b[a]-d[c])/b[a])>0.01))+"\t")