\chapter{Page Rank}
In questo capitolo ci porremo il problema di calcolare il \emph{page rank} di un determinato
insieme di pagine Web (o, analogamente, documenti \emph{ipertestuali}). \\
Con Page Rank intendiamo un ``valutazione'' data ad una determinata pagina o documento con lo scopo
di stimare la sua importanza\footnote{che è ovviamente un concetto poco matematico, ma ci occuperemo
di formalizzarlo meglio in seguito.}. Un algoritmo derivato da quello che verrà in seguito
esposto è utilizzato attualmente da Google (\href{http://www.google.it}{http://www.google.it})
per ordinare i risultati di una query nel suo motore di ricerca, così come da tutti i principali
\emph{search engine} presenti sul Web.
 
\section{Presentazione del modello}
Prima di tutto dobbiamo occuparci di definire il modello matematico che vogliamo usare
per schematizzare la nostra valutazione di importanza. \\
In generale l'idea che sta alla base del Page Rank è dare importanza ad una pagina
basandosi sui link\footnote{Un link in una pagina Web, ma anche in un documento PDF come questo,
è un ``puntatore'' ad un altro documento. In genere l'utente può attivare il puntatore semplicemente
cliccandoci sopra con il mouse} che in essa sono presenti.
 
Alcune idee che potrebbero sembrare valide potrebbero essere valutare l'importanza di una pagina
in base al numero di link presenti nella pagina, oppure al numero di link che puntano ad essa.
 
Dopo una breve riflessione ci si rende conto che entrambi i metodi introducono svariate problematiche
di affidabilità. In generale, infatti, il creatore di una pagina ha il controllo di quanti link
sono presenti in una pagina, e non fatica molto a creare delle pagine fittizie che puntano a quella
appena creata. La valutazione sarebbe quindi falsata e molto vulnerabile a frode.
 
Si può però cercare di affinare l'idea di base in modo da risolvere questi inconvenienti.
Per convenzione assumiamo di numerare tutte le pagine presenti nel Web con dei numeri naturali
$1 \ldots N$ (dove $N$ è un numero dell'ordine di $10^{10}$).
Fissata una certa pagina $j$ consideriamo i link che puntano a lei. Risulta chiaro che non solo
il \textbf{numero} dei link è importante per valutare la sua importanza, ma a sua volta anche l'\textbf{importanza}
di questi dovrebbe influire. Per esempio, se la homepage di Yahoo! punta
alla pagina di una data persona, è naturale che questa pagina acquisisca una certa rilevanza. Viceversa, se sulla
mia pagina personale metto un link al sito della stessa persona l'importanza di questo varia di poco\footnote{
a meno che nel momento in cui voi leggiate queste pagine io non sia diventato una persona con una certa influenza,
ma questo è probabilmente un caso da scartare}.
Un'ultima e non meno importante considerazione è che se anche Yahoo! punta al mio sito ma lo fa in una pagina
con un altro centinaio di link è una faccenda molto diversa dall'essere puntati dall'unico link in homepage.
 
Stiamo quindi arrivando ad un altro modello, in cui le pagine \emph{distribuiscono la loro importanza fra i link al
loro interno}. Questo modello è molto più ragionevole del precedente, ed è molto meno predisposto a dare
risultati falsati. Vorremo però averne una descrizione matematica precisa.
 
\section{Un po' di matematica}
Consideriamo la matrice quadrata $H = (h_{ij})$, che chiameremo \emph{matrice delle connessioni}, definita
nel seguente modo
\[
 h_{ij} = \left\{ \begin{array}{ll}
                   1 & \text{se esiste un link dalla pagina} \ i \ \text{alla pagina} \ j \\
		   0 & \text{altrimenti}
                  \end{array} \right.
\]
Questa matrice avrà dimensione $N \times N$ ma sarà sparsa, ovvero la maggior parte dei suoi valori saranno $0$.
Per fissare le idee consideriamo la Figura~\ref{fig:web}, un modello (molto) semplificato del Web, in cui esistono solamente
3 pagine.
\begin{figure}[hb]
\[
 \xymatrix{
{\bullet{1}} \ar@/_/[dr] \ar@/^1pc/[drr] & \\
 & {\bullet{2}} \ar@/_/[ul]& {\bullet{3}} \ar@/^/[l]}
\]
\caption{Esempio di una rete composta da 3 pagine}
\label{fig:web}
\end{figure}
 
La matrice di connessione $H$ associata a questo grafo è la seguente
\[
 H = \left[ \begin{array}{ccc}
             0 & 1 & 1 \\
	     1 & 0 & 0 \\
	     0 & 1 & 0
            \end{array} \right]
\]
\begin{os} \label{os:pageranknumcoll}
 \`E facile verificare che valgono le seguenti uguaglianze
 \[
 \left\{ \begin{array}{ll}
  \sharp \{ \text{collegamenti uscenti da} \ i \} & = \sum_{j=1}^{N} h_{ij} \\
  \sharp \{ \text{collegamente che arrivano a} \ j \} & = \sum_{i=1}^{N} h_{ij}
  \end{array} \right.
 \]
\end{os}
In particolare sarebbe piuttosto immediato calcolare, partendo dalla matrice $H$, i valori
del Page Rank con i metodi ipotizzati all'inizio del capitolo (per altro reputati incompleti)
 
Siamo dunque interessati a calcolare il vettore $W$ dove $w_i$ è l'importanza della pagina $i$,
in modo che rispetti le linee guide definite precedentemente.
 
Chiamiamo $d_i$ il numero di link uscenti dalla pagina $i$
(facilmente calcolabile come visto nell'Osservazione~\ref{os:pageranknumcoll}). Chiameremo inoltre
$D$ la matrice con i $d_i$ sulla diagonale e nulla per i restanti valori.
Introdotta questa notazione, secondo il nostro modello ogni componente di $W$ dovrebbe rispettare
la seguente equazione
\[
 w_j = \sum_{i=1}^{N} w_i \frac{h_{ij}}{d_i}
\]
che rappresenta esattamente l'idea che l'importanza di ogni pagina è la somma delle importanze di quelle
che la linkano, in cui ogni addendo è normalizzato sul numero di link uscenti\footnote{probabilmente l'equazione
è molto più chiara di questa frase, ma non commentare mi sarebbe sembrato inopportuno}.
\`E quindi chiaro come il nostro problema si sia trasformato in un problema di autovalori ed autovettori. In
particolare, se poniamo $M = D^{-1}H$ la nostra richiesta si può scrivere nella forma più compatta
$\trasp{W} = \trasp{W}M$, ovvero $\trasp{W}$ è un autovettore sinistro di $M$. Sorgono spontanee alcune
domande a cui dobbiamo rispondere prima di procedere oltre.
\begin{enumerate}[(a)]
 \item \label{en:pagerank:a} Cosa succede se per qualche $i$ si ha $d_i = 0$? Stiamo tacitamente assumendo che $D$ sia invertibile
 ma questo non è garantito a priori;
 \item \label{en:pagerank:b} L'autovalore $1$ fa parte dello spettro delle matrice $M$? Se questo non fosse vero la nostra
 richiesta non ammette soluzione, e quindi possiamo rinunciare in partenza alla nostra idea;
 \item \label{en:pagerank:c} Se l'autovettore $W$ esiste, è unico\footnote{In realtà ci stiamo chiedendo se il suo
autospazio ha dimensione $1$. Siamo ben consci che se ne troviamo uno ogni suo multiplo scalare continuerà
ad essere un autovettore. Per altro, moltiplicare le valutazioni omogeneamente per lo stesso cofficiente non
sarebbe in nessun modo un problema}?. Se così non fosse
 sorgerebbero dei seri problemi su come interpretare i risultati ottenuti;
 \item \label{en:pagerank:d} Possiamo scegliere l'autovettore $W$ in modo che $w_i \geq 0$ per ogni $i$? Questo non è necessario,
 a priori, ma sembrerebbe coerente fare in modo che delle votazioni siano numeri positivi;
\end{enumerate}
 
Cominciamo col rispondere alla domanda~(\ref{en:pagerank:a}). Ovviamente potrebbe benissimo succedere
che per qualche $i$ $d_i = 0$. Questi ``nodi'' vengono chiamati \emph{dangling node}, ovvero ``nodi penzolanti'',
che non portano a niente. Per risolvere questo problema si sostituisce la matrice $H$ con una matrice $\hat H$
tale che se $h_{ij} = 0 \: \forall j$ allora $\hat h_{ij} = 1 \: \forall j$, e in caso contrario $\hat h_{ij} = h_{ij}$.
Questa matrice risolve evidentemente il problema e cambia di poco il modello. Supporre che una determinata pagina
non distribuisca la sua importanza a nessuno o la distribuisca in maniera analoga a tutti non è un gran cambiamento.
 
La domanda~(\ref{en:pagerank:b}) ha invece una risposta immediata. Sì, l'autovalore $1$ fa parte dello spettro
ed un autovettore (destro) a lui associato è quello composto di tutti $1$. La verifica è un banale conto.
 
Per rispondere alle domande~(\ref{en:pagerank:c}) e (\ref{en:pagerank:d}) ci appoggiamo ad un risultato
che non dimostreremo (in questa sede) ottenuto da Perron-Frobenius sulla Teoria delle Matrici non negative.
\begin{te}[di Perron-Frobenius] \label{te:PerronFrobenius}
 Se $A$ è irriducibile allora $\rs{A} > 0$ è un autovalore semplice ed
 esistono $x,y$ autovettori destro e sinistro
 con componenti strettamente positive.
\end{te}
\begin{os}
 Sembrerebbe, a priori, che non siamo in grado di determinare il raggio spettrale di $A$. Consideriamo però
 che $1$ è nello spettro e che $||A||_{\infty} = 1$ e ogni norma indotta supera il raggio spettrale. Si può
 quindi concludere che $\rs{A} = 1$ e applicare il Teorema~\ref{te:PerronFrobenius}.
\end{os}
Sfortunatamente si può osservare che, pur con tutte queste buone proprietà, in generale i metodi applicati
per trovare l'autovettore (il metodo delle potenze) potrebbe non essere convergente, perché gli autovalori
diversi da $1$ potrebbero ancora avere modulo $1$\footnote{Come visto nella teoria, la velocità di convergenza
del metodo di iterazione è data dal rapporto fra il più grande autovalore diverso da quello di cui si sta cercando
l'autovettore e l'autovalore in questione. In questo caso il rapporto potrebbe essere $1$, che darebbe una situazione
di non convergenza}.
Introduciamo una definizione che ci permetterà di formalizzare meglio il problema.
\begin{de}
  Una matrice $A$ si dice \emph{ciclica} se esiste un sottoinsieme dei suoi autovalori $J \subseteq \spe{A}$
  tale che $\forall \lambda,\mu \in J$ si ha che $|\lambda|=|\mu|$.
\end{de}
Per superare questa difficoltà si può, ancora una volta, modificare il modello. Fissiamo $\gamma > 0$ e consideriamo
che il visitatore
possa in ogni momento, con probabilità $1-\gamma$, decidere di aprire una pagina a caso. Potremmo rappresentare
questa nuova situazione sostituendo $M$ con la matrice
\[
 A = \gamma M + (1-\gamma) \left[ \begin{array}{c}
                                   1 \\
				   \vdots \\
				   1
                                  \end{array} \right] \left[ 1 \: \hdots \:  1 \right]
\]
ovvero con una combinazione convessa della matrice $M$ con la matrice di soli $1$. Il valore di $\gamma$ si può
scegliere a piacimento; Google usa un valore $\gamma = 0.85$. In questo caso la matrice $A$ non ha elementi nulli
e quindi è irriducibile e non è ciclica\footnote{La condizione di ciclicità è quella che crea problemi
con gli autovalori di modulo uguale al raggio spettrale}. Possiamo infine utilizzare un'altra versione del Teorema~\ref{te:PerronFrobenius}
che dice
\begin{te}
 Sia $A$ una matrice con elementi strettamente positivi. Allora per ogni suo autovalore $\lambda \neq \rs{A}$ si ha
 che vale la diseguaglianza $|\lambda| < |\rs{A}|$.
\end{te}
Quest'ultimo ci assicura la convergenza del nostro metodo.
 
\subsection{Implementazione del metodo delle potenze}
Presenteremo in questa sezione una versione del metodo delle potenze adattata alle nostre esigenze.
Consideriamo la matrice $A$ di cui vogliamo trovare l'autovettore sinistro relativo all'autovalore $1$.
Sappiamo che esiste la forma di Jordan di $A$ ed in particolare una matrice $S$ tale che
\[
 S^{-1} A S = J \qquad J = \left[ \begin{array}{c|cc}
                                   1 & \multicolumn{2}{|c}{\herm{0}} \\ \hline
               \multirow{2}{*}{$0$} & \sblocke{\hat J}{2} \\
                                   & &
                                  \end{array} \right]
\]
Osserviamo ora che $e_1$ è autovettore destro di $J$ ed anche autovettore sinistro. Entrambi sono relativi
all'autovalore $1$. Sappiamo quindi che $e = Se_1$ è autovalore destro di $A$ e quindi è il vettore
con tutte le componenti uguali ad $1$, mentre $\trasp{w} = \trasp{e_1}S^{-1}$ è autovettore sinistro di $A$,
ovvero il vettore che stiamo cercando.
In particolare quanto visto nella precedente sezione ci assicura che $\rs{\hat J} < 1$ e quindi possiamo osservare
che
\[
 \lim_{k\to\infty} J^{k} = e_1\trasp{e_1} \quad \text{e quindi} \quad \lim_{k\to\infty} A^{k} =
 \lim_{k\to\infty} SJ^kS^{-1} = Se_1 \trasp{e_1}S^{-1} = e\trasp{w}
\]
Questo (apparentemente) risolve il nostro problema in quanto $e\trasp{w}$ non è altro che la matrice le cui
righe coincidono con il vettore $\trasp{w}$ che è precisamente quanto stiamo cercando. \\
La complessità di questo calcolo non è però trascurabile, in quanto la moltiplicazione di una matrice
costa in generale $O(n^3)$, decisamente troppo se $n \approx 10^{10}$ come nel nostro caso.
\begin{os}
 Consideriamo dunque un vettore $x_0$ tale che $\trasp{e}x_0 = 1 = ||x_0||_{1}$, e definiamo la successione per ricorrenza
 \[
  \left\{ \begin{array}{ll}
           \trasp{x_0} & = \trasp{x_0} \\
           \trasp{x_{k+1}} & = \trasp{x_{k}} A
          \end{array} \right.
 \]
Possiamo notare che in generale $x_k = x_0 A^{k}$ e quindi data la nostra scelta oculata di $x_0$ abbiamo
che
\[
 \lim_{k\to\infty} \trasp{x_k} = \lim_{k\to\infty} \trasp{x_0} A^{k} = x_0 e \trasp{w} = \trasp{w}
\]
In questo modo siamo scesi da una complessità $O(n^3)$ ad una $O(n^2)$.
\end{os}
Sfortunatamente $O(n^2)$ nel nostro caso continua ad essere un numero troppo grande; arrivati a questo
punto ci accorgiamo che possiamo sfruttare la particolarità di $A$, che non è assolutamente
una matrice qualsiasi! Sappiamo infatti che è composta quasi esclusivamente di zeri, e quindi una
memorizzazione efficiente ci permetterebbe anche un'implementazione efficiente dei procedimenti
di moltiplicazione.
%% TODO: QUalcosa di più sull'implementazione del metodo
 
\subsection{Note tecniche}
Analizzato il problema matematicamente, sorge il bisogno di analizzare le difficoltà di implementazione dovute
all'hardware attualmente disponibile. Una matrice delle dimensioni di $A$ non è realmente memorizzabile
elemento per elemento su nessun calcolatore, e non è pensabile nemmeno in alcun tipo di cluster. \\
Un rapido conto mostra come, anche ipotizzando di usare dei tipi \verb-float- a 4 byte (in ogni
caso presumibilmente troppo imprecisi) otteniamo uno spazio necessario di $N^2 \cdot 4$ byte, ovvero
circa $10^{20}$ byte = $10^{11}$ GB di RAM, probabilmente un quantitativo superiore a quella di tutti i computer
del mondo messi assieme\footnote{può darsi che nei prossimi anni quest'affermazione smetta di essere vera, ma
non penso sia rilevante ai fini del problema}.
Dobbiamo quindi adottare degli stratagemmi per gestire la memoria in modo efficiente.
 
Una prima osservazione è che la nostra matrice $A$ si può scrivere nel seguente modo
\[
A = \gamma D^{-1} \hat H + (1-\gamma) e\trasp{v}
\]
dove $\gamma$ e $w$ sono rispettivamente il parametro ed il vettore di personalizzazione. Per memorizzare
$A$ possiamo quindi memorizzare $D$ ed $\hat H$, $\gamma$ e $w$. \\
Questo ci dà dei vantaggi perché la matrice $\hat H$ è ``quasi nulla'' nel senso che è fatta solamente di $0$ ed $1$,
ma la percentuale di $0$ è molto alta\footnote{è infatti presumibile che il numero di link uscenti da una pagina
sia trascurabile rispetto al numero di pagine esistenti}. Possiamo quindi limitarci a memorizzare i seguenti dati: \\[5pt]
\begin{tabular}{|lll|} \hline
 \textbf{Variabile} & \textbf{Descrizione} & \textbf{Spazio occupato} \\ \hline \hline
 $N$ & l'ordine della matrice $A$, ovvero il numero di pagine esistenti & 4 byte \\
 $m$ & il numero di campi di $\hat H$ diversi da $0$ & \textasciitilde $N$ byte \\
 $i, j$ & vettori delle coordinate dei campi diversi da $0$ di $\hat H$ & \textasciitilde $N$ byte \\
 $\gamma$ & coefficiente di personalizzazione & 4 byte \\
 $w$ & Vettore di personalizzazione & \textasciitilde $N$ byte \\
 $d$ & il vettore con gli elementi diagonali di $D$ & \textasciitilde $N$ byte \\
\hline \end{tabular} \\
 
Possiamo ora porci il problema di calcolare $d$. Questo è reso apparentemente più complicato dal modo in cui abbiamo salvato
la matrice. In realtà possiamo usare un algoritmo di questo tipo
 
\begin{lstlisting}[frame=tb,caption=Calcolo di $d$,label=lst:calcolovettorepers]{}
 ! Inizializziamo il vettore d a 0
 d = 0;
 ! Questo ciclo calcola i link uscenti da i per ogni i
 DO k=1,m
   d(i(k)) = d(i(k)) + 1
 END DO
\end{lstlisting}
 
\begin{os}
 L'aver scelto questo tipo di memorizzazione per la matrice non solo non complica i generali procedimenti
 di moltiplicazione, ma riduce addirittura la complessità di questi perché iteriamo sugli elementi diversi da
 $0$ invece che su tutti.
\end{os}
 
 
\section{Una prima implementazione}
In questa sezione ci occuperemo di presentare una prima implementazione della soluzione al problema.
Non ci preoccuperemo troppo di ottimizzare il codice, compito che affronteremo in seguito.
Esporremo inizialmente la struttura generale del programma e poi costruiremo tutte le subroutine necessarie.
 
\subsection{La struttura}
Il programma dovrà comporsi delle seguenti parti
\begin{description}
 \item[Lettura della matrice] Supporremo di avere un file da cui leggere i dati della matrice. Dobbiamo
 definire la sintassi del file e scrivere un subroutine per caricarlo in memoria;
 \item[Calcolo del vettore d] Data la matrice di connessione $H$ vorremo calcolare la somma di ogni
 riga che ci servirà in seguito per rinormalizzare il vettore di iterazione;
 \item[Iterazione] Generiamo un vettore $x$ casuale e lo riscaliamo in modo che $\trasp{e} x = 1$. Fatto
 questo cominciamo a calcolare l'iterazione, fino a che l'errore non è piccolo o ragguingiamo il massimo
 numero di iterazioni impostato;
 \item[Salvataggio] Salviamo su un file il vettore calcolato;
\end{description}
 
Per questo avremo bisogno anche di alcune subroutine, ad esempio quella che ci calcoli il prodotto matrice
vettore data la matrice sparsa come l'abbiamo intesa noi.
 
\subsection{Prodotto di matrice sparsa per un vettore}
Esponiamo qui un'idea di codice per calcolare questo prodotto matrice vettore. Ricordiamo che per noi una
matrice sparsa è composta di due scalari \lstinline-n-,\lstinline-m- (rispettivamente dimensione e numero di elementi diversi da $0$)
e due vettori \lstinline-i-, \lstinline-j- di dimensione $m$ che sono le coordinate dei punti che valgono $1$.\\
Il codice è il seguente
\begin{lstlisting}[frame=tb,caption=Calcolo del prodotto matrice vettore,label=lst:matrixvec]
 SUBROUTINE dotprod(n,m,i,j,d,x,v,y)
  IMPLICIT NONE
  ! Il parametro dp dei tipi a doppia precisione
  INTEGER, PARAMETER :: dp = KIND(0.D0)
 
  ! Dichiarazione degli argomenti della funzione
  INTEGER, INTENT(IN) :: n,m
  INTEGER, DIMENSION(m), INTENT(IN) :: i,j
  INTEGER, DIMENSION(n), INTENT(IN) :: d
  REAL(dp), DIMENSION(n), INTENT(IN) :: x, v
 
  ! Dato di output
  REAL(dp), DIMENSION(n) :: y
 
  ! Alcune variabili di comodo
  REAL(dp), DIMENSION(n) :: xx
  REAL(dp) :: s, gamma = 0.85
  INTEGER :: k
 
  ! Calcoliamo x trasposto per la matrice D con i d(i) sulla
  ! diagonale,  e salviamolo in xx, in modo da non modificare x
  ! Ricordiamo che i dati in fortrano vengono sempre passati per
  ! riferimento! Se modificassimo x qui dentro risulterebbe cambiato
  ! anche all'esterno.
  DO k=1,m
   IF(d(k)/=0) THEN
     xx(i(k)) = x(i(k))/d(k)
   ELSE
     xx(i(k)) = x(i(k))/n
   END IF
  END DO
 
  ! Ora eseguiamo la moltiplicazione di xx per la parte non dangling
  ! della matrice H cappuccio e salviamo il risultato in y
  y = 0.D0
  DO k=1,m
    y(j(k)) = y(j(k)) + xx(i(k))
  END DO
 
  ! Ora calcoliamo le righe dangling dalle quali prima non siamo "passati"
  ! perche' i(k) scorre solo sulle righe dove ci sono elementi non nulli
  DO k=1,n
    IF(d(k) == 0)
      s = s + xx(k)
    END IF
  END DO
 
  ! Sommiamo i due contributi
  y = y + s
 
  ! Calcolo la personalizzazione
  y = y * gamma + (1-gamma)/n
 
 END SUBROUTINE dotprod
\end{lstlisting}
 
%TODO: Presentare un'implementazione completa del programma