Page Rank

Let’s use the 6 page universe that we had in the course notes. For this directed graph, perform the following calculations in R.

Form the A matrix. Then, introduce decay and form the B matrix as we did in the course notes

\[ A = \begin{bmatrix} 0 & 0 & \frac{1}{4} & 0 & 0 & 0 \\ \frac{1}{2} & 0 & \frac{1}{4} & 0 & 0 & 0 \\ \frac{1}{2} & 1 & 0 & 0 & 0 & \frac{1}{2} \\ 0 & 0 & 0 & 0 & \frac{1}{2} & \frac{1}{2} \\ 0 & 0 & \frac{1}{4} & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{4} & \frac{1}{2} & \frac{1}{2} & 0 \end{bmatrix} \]

(A <- matrix (c(  0, 0,  1/4,  0,   0,   0,
               1/2, 0,  1/4,  0,   0,   0,
               1/2, 1,    0,  0,   0, 1/2,
                 0, 0,    0,  0, 1/2, 1/2,
                 0, 0,  1/4, 1/2,  0,   0,
                 0, 0,  1/4, 1/2,1/2,   0), byrow=TRUE, nrow=6))
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  0.0    0 0.25  0.0  0.0  0.0
## [2,]  0.5    0 0.25  0.0  0.0  0.0
## [3,]  0.5    1 0.00  0.0  0.0  0.5
## [4,]  0.0    0 0.00  0.0  0.5  0.5
## [5,]  0.0    0 0.25  0.5  0.0  0.0
## [6,]  0.0    0 0.25  0.5  0.5  0.0

Introduce decay and form the B matrix

d <- 0.85
nrows <- nrow(A)
(B <- d * A + ( (1-d) / nrows ))
##       [,1]  [,2]   [,3]  [,4]  [,5]  [,6]
## [1,] 0.025 0.025 0.2375 0.025 0.025 0.025
## [2,] 0.450 0.025 0.2375 0.025 0.025 0.025
## [3,] 0.450 0.875 0.0250 0.025 0.025 0.450
## [4,] 0.025 0.025 0.0250 0.025 0.450 0.450
## [5,] 0.025 0.025 0.2375 0.450 0.025 0.025
## [6,] 0.025 0.025 0.2375 0.450 0.450 0.025

Start with a uniform rank vector r and perform power iterations on B till convergence.

#The uniform rank vector for the above 6 page universe [i.e., each carrying the equal probability ]
(r <- rep(1/nrows, nrows))
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

Power Iterations

Power Iteration Formula: \[ r_f = A^n \times r_i \] (here, \(n\) is the number of iterations of the algorithm that we want to run. As we run through these iterations, the ranks of the pages converge and stabilize).

powerIteration <- function(A, r, n) {
  #This function performs the power iterations to calculate the probability of landing on a given page.
  #Inputs:
  # A = Input matrix
  # r = initial rank of the urls
  # n = Number of Iterations
  #
  #Output:
  #
  # Final probability matrix computed

  Iter = diag(dim(A)[1])
  for ( i in 1:n)
  {
    Iter = Iter %*% A
  }
  
  return (Iter %*% r)
  
}

Now, lets call the above function with B, r for few timess

(powerIteration(B, r, 20))
##            [,1]
## [1,] 0.07735889
## [2,] 0.11023642
## [3,] 0.24639470
## [4,] 0.18635383
## [5,] 0.15655925
## [6,] 0.22309691
(powerIteration(B, r, 30))
##            [,1]
## [1,] 0.07735886
## [2,] 0.11023638
## [3,] 0.24639464
## [4,] 0.18635389
## [5,] 0.15655927
## [6,] 0.22309696
(r_j <- powerIteration(B, r, 40))
##            [,1]
## [1,] 0.07735886
## [2,] 0.11023638
## [3,] 0.24639464
## [4,] 0.18635389
## [5,] 0.15655927
## [6,] 0.22309696

Notice that the convergence occurs after 30 iterations !

Eigen-decomposition of \(B\)

Compute the eigen-decomposition of B and verify that you indeed get an eigenvalue of 1 as the largest eigenvalue and that its corresponding eigenvector is the same vector that you obtained in the previous power iteration method. Further, this eigenvector has all positive entries and it sums to 1.

eigenDecom <- eigen(B)
(maxEigenVal <- max(as.numeric(eigenDecom$values)))
## Warning: imaginary parts discarded in coercion
## [1] 1
#Find the corresponding eigen vector
ndindexes <- which(round(maxEigenVal, 2) == 1.00)
matrix(eigenDecom$vectors[,ndindexes], nrow=6, byrow=TRUE)
##               [,1]
## [1,] -0.1784825+0i
## [2,] -0.2543376+0i
## [3,] -0.5684822+0i
## [4,] -0.4299561+0i
## [5,] -0.3612138+0i
## [6,] -0.5147297+0i

This is not same as the one we exepcted in the power iteration method (r_j)! , eigen function must be normalizing the eigenvector. Let’s try to find the unit vector for (r_j)

unitVector <- function(v) {
  vlen <- sqrt(sum(v^2))
  return ( v / vlen)
}

(unitVector(r_j))
##           [,1]
## [1,] 0.1784825
## [2,] 0.2543376
## [3,] 0.5684822
## [4,] 0.4299561
## [5,] 0.3612138
## [6,] 0.5147297

Use the igraph package in R and its page.rank method to compute the Page Rank of the graph as given in A.

library(igraph)
## 
## Attaching package: 'igraph'
## 
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## 
## The following object is masked from 'package:base':
## 
##     union
#Create adj graph from the given matrix A.  IMPORTANT NOTE: Since we set each column as vertex, we need to transpose the matrix A.
graphObj <- graph.adjacency(t(A), weighted = TRUE, mode = "directed")
plot(graphObj)

Let’s compute the page rank using the igraph

(prVec <- page.rank(graphObj)$vector)
## [1] 0.07735886 0.11023638 0.24639464 0.18635389 0.15655927 0.22309696

Let’s compute the unit vector to compare with our earlier results

(matrix(unitVector(prVec), nrow=6, byrow=TRUE))
##           [,1]
## [1,] 0.1784825
## [2,] 0.2543376
## [3,] 0.5684822
## [4,] 0.4299561
## [5,] 0.3612138
## [6,] 0.5147297
round(prVec) == round(r_j)
##      [,1]
## [1,] TRUE
## [2,] TRUE
## [3,] TRUE
## [4,] TRUE
## [5,] TRUE
## [6,] TRUE