Graphs, adjacency matrices, and pagerank.

Introduction

In this notebook, I discuss four approaches to pagerank.  I begin by illustrating some of the tools built in to Mathematica to handle graphs.  I use these tools to write a Monte Carlo simulation of a websurfer and produce the pageranks.  That’s the first approach.

Then, do a direct computation by solving the system of equations satsified by pageranks.  That’s the second approach.

Then, I use some tools to extract a matrix, I call it c, from the graph.  The vector of pageranks is an eigenvector of the matrix c with eigenvalue 1.  I then find that eigenvector in two ways.  The first is with a built in Mathematica command (the third approach) and the second way is with a simple numerical method called the “Power Method” (the fourth approach).

Graphs

Mathematica has a lot of graph theoretic capabilities built in.  Mathematica has features to deal with many kinds of graphs (undirected, directed, multigraphs, ...) but I’ll focus on directed graphs since that’s what’s relevant for right now.  

For example, here’s a way to input a directed graph that you can think of as a toy model of linked webpages:

g=Graph[{
1\[DirectedEdge]2,3\[DirectedEdge]2,4->1,5->6,6->1,1\[DirectedEdge]7,
2\[DirectedEdge]6,2\[DirectedEdge]7,
3->5,4->3,4->5,4->8,
5->8,
6->8,
7->4,7->5,7->6,
8->2,8->7
}]

Here’s a labelled picture of our graph g:

GraphPlot[g,VertexLabeling->True]

Four ways to find the page ranks of a directed graph

A Monte Carlo simulation

We’ll write a program that inputs a graph and a number of iterations and outputs and outputs the proportion of time a random websurfer spends on each vertex of the graph:

First, let’s see how we can extract the information we need from the graph.  There are simple commands to get the list of vertices and the list of edges.

VertexList[g]

EdgeList[g]

If we want particular edges that match a pattern, say for example, we want all edges that originate with vertex 1, we’d use:

EdgeList[g,1\[DirectedEdge]_Integer]

If we want the “second” edge coming out of vertex 1, we’d use:

EdgeList[g,1\[DirectedEdge]_Integer][[2]]

And if we want the vertex that the second edge coming out of vertex 1 terminates with, we’d use:

EdgeList[g,1\[DirectedEdge]_Integer][[2,2]]

Here’s a random choice of a vertex that vertex 4 connects to:

RandomChoice[EdgeList[g,4\[DirectedEdge]_Integer]][[2]]

websim[graph_,iterations_]:=

Module[
{vertices,n,vertexcount,currentvertex,newvertex}, (* here are the variables we will use in our program *)

vertices=VertexList[graph]; (* get the vertices of the graph *) 
n=Length[VertexList[graph]];  (* set n= the number of vertices  *)
 
vertexcount=Table[0,{i,1,n}]; (* set the number of visits to each vertex equal to zero *)

currentvertex=RandomChoice[vertices];  (* choose an initial vertex randomly *)

Do[
newvertex=
RandomChoice[EdgeList[graph,currentvertex\[DirectedEdge]_Integer]][[2]];  (*choose a new vertex randomly among the ones that the current vertex is connected to *)

currentvertex=newvertex;  
vertexcount[[currentvertex]]++; ,(* make the new one the current one and increment the number of visits to the current page *)

{iterations}  (* Loop "iterations" number of times *)
];

vertexcount/iterations  (* output the proportion of visits to each webpage *)

]

Let’s simulate a surfer visiting one hundred thousand pages:

pagerankdata=websim[g,100000];

Let’s see the results:

TableForm[N[pagerankdata]]

Eigenvalue method

The Adjacency Matrix

Mathematica can compute the a adjacency matrix of a graph.  Since adjacency matrices for graphs can be quite large, and because they have many entries that are zero, Mathematica will treat an adjacency matrix as a “Sparse Array” by default.  That means that that it’s only storing the nonzero entries (and not wasting time or space storing zeros).

AdjacencyMatrix[g]

However, we can easily convert the sparse array to a normal matrix (which in Mathematica is a list of lists, which are the rows of the given matrix).

a=Normal[AdjacencyMatrix[g]]

And if we want, we can look at the matrix arranged in rows and columns:

MatrixForm[a]

The i-th row of this matrix is the i-th list in a.  For example, the fourth row is the fourth list in the list a:

a[[4]]

The entry in the 4th row and the 3rd column of the matrix a is the 3rd element of the 4th row.  There are two ways to get this element from a:

a[[4]][[3]]

a[[4,3]]

It will be helpful to calrify the convention for the adjacency matrix:  There is a 1 in row i, column j if there is a directed edge starting at vertex i and terminating at vertex j.  So, the entry “1” in the second row, sixth column is there because there is an edge from vertex 2 to vertex 6.

For the pagerank algorithm, we’d like a matrix with the entry 1/(Subscript[d, j]) in the i-th row and j-th column if there is an edge from vertex j to vertex i, and vertex j has d outgoing edges. 

We can determine the number of outgoing edges for each vertex:  The number of edges coming out of vertex i is the total number of 1’s in row i.  For example, the number edges originating at vertex 7 is

Total[a[[7]]]

Here is a list of degrees:

d=Table[Total[a[[i]]],{i,1,8}]

So, here is the matrix of interest for the pagerank.

c=Table[a[[j,i]]/d[[j]],{i,1,8},{j,1,8}];
MatrixForm[c]

Direct approach to eigenvalue problem

Note that the equation c.probabilities==probabilities says that the vector of probabilities  is an eigenvector of the matrix c with eigenvalue 1.  Mathematica has built in functions to compute eigenvectors, let’s use them:

First, as an 8x8 matrix, c may have as many as eight distinct eigenvectors corresponding to eight distinct eigenvalues.  This is in fact the case (some are complex)

The list of eigenvectors is long, and most are expressed in terms of the solutions to various polynomial equations, which are not solved exactly.

Eigenvectors[c]

N[Eigenvalues[c]]

The one we are interested in is the first :

Eigenvectors[c][[1]]

Any multiple of this vector is also an eigenvector.  The multiple we’re interested in, the one that gives the pagerank probabilities, is the multiple for which the sum of the entries equals 1.

eigenvectorsolution=Eigenvectors[c][[1]]/Total[Eigenvectors[c][[1]]]

TableForm[100N[eigenvectorsolution]]

Power method approach to eigenvalue problem

We will go into the theory behind the “power method” in class, but by way of introduction, here is a brief illustration:

We will choose a vector randomly, and iteratately apply the matrix c to that vector.  The result approaches an eigenvector that corresponds to the largest eigenvalue.

multiplybyc[vector_]:=c.vector

myrandomvector=RandomReal[1,8]

c.c.c.myrandomvector

Nest[multiplybyc,myrandomvector,3]

Nest[multiplybyc,myrandomvector,20]

approxeigenvector=%/Total[%]

Direct computation

Let (Subscript[p, i])be the probability of being at vertex i.  We compute the probability (Subscript[p, i])satisfies 

(Subscript[p, i])=∑(1/(Subscript[d, j]))(Subscript[p, j])

Where the sum is over vertices j with an edge that terminates in vertex i.  

This gives the following system of equations:

equations=Table[(Subscript[p, i])==Sum[a[[j,i]]/d[[j]](Subscript[p, j]),{j,1,8}],{i,1,8}];
TableForm[equations]

This system has infinitely many solutions, but if we add the equation (Subscript[p, 1])+(Subscript[p, 2])+(Subscript[p, 3])+(Subscript[p, 4])+(Subscript[p, 5])+(Subscript[p, 6])+(Subscript[p, 7])+(Subscript[p, 8])==1 we get a system of equations with a unique solution, which we can solve:

equations=Append[equations,Sum[(Subscript[p, i]),{i,1,8}]==1];
TableForm[equations]

probabilities={(Subscript[p, 1]),(Subscript[p, 2]),(Subscript[p, 3]),(Subscript[p, 4]),(Subscript[p, 5]),(Subscript[p, 6]),(Subscript[p, 7]),(Subscript[p, 8])}

solution=Solve[equations,probabilities]

TableForm[ Table[{i,100 N[solution[[1,i,2]]]},{i,1,8}]]

Notice that the direct computation of the probabilites involved solving a system of equations.  This system of equaitons (without the normalization (Subscript[p, 1])+(Subscript[p, 2])+(Subscript[p, 3])+(Subscript[p, 4])+(Subscript[p, 5])+(Subscript[p, 6])+(Subscript[p, 7])+(Subscript[p, 8])==1 ) is equivalent to the single matrix equation where c is the matrix defined earlier

c.probabilities==probabilities

Comparison of Results

Let’s pause to compare our four soluitons to the pagerank problem:

table=N[ Table[{i,100 pagerankdata[[i]],100 solution[[1,i,2]],100N[eigenvectorsolution[[i]]],100 approxeigenvector[[i]]},{i,1,8}]];

TableForm[Prepend[table,{"vertex", "monte carlo", "direct solution", "direct eigenvalue","approx eigenvalue"}]]

Timing[websim[g,1000000];]

Timing[Solve[equations,probabilities]]

Timing[Eigenvectors[c][[1]]]

Timing[Nest[multiplybyc,myrandomvector,20]]

Exercises

Problem 1

Try to construct a graph with 3 vertices so that the pagerank so that the pageranks of the vertices are (1/3, 1/3, 1/3).  

Try to construct a graph with 3 vertices so that the pagerank so that the pageranks of the vertices are (1/6, 1/3, 1/2).

Problem 2

As a test of your understanding of how Mathematica treats graphs, write your own Monte Carlo program that inputs a graph g and a number of iterations, and outputs the pageranks.

As a test of your understanding of how Mathematica treats matrices, write your own Monte Carlo program that inputs a matrix that is the adjacency matrix of a directed graph and a number of iterations, and outputs the pageranks of the graph whose adjacency matrix was given.