PLAYING WITH PAGERANK

(1) Form the A matrix. Then introduce decay and form the B matrix as in course notes.

options(digits=2)
# Matrix A formed
A <- matrix(c(0, 1/2, 1/2, 0, 0, 0,
              0, 0, 0, 0, 0, 0,
              1/3, 1/3, 0, 0, 1/3, 0,
              0, 0, 0, 0, 1/2, 1/2,
              0,0, 0, 1/2, 0, 1/2,
              0, 0, 0, 1, 0, 0), nrow=6)


print(A)
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  0.0    0 0.33  0.0  0.0    0
## [2,]  0.5    0 0.33  0.0  0.0    0
## [3,]  0.5    0 0.00  0.0  0.0    0
## [4,]  0.0    0 0.00  0.0  0.5    1
## [5,]  0.0    0 0.33  0.5  0.0    0
## [6,]  0.0    0 0.00  0.5  0.5    0
# Decay value applied and matrix B is formed
d <- 0.85
B <- d * A + (0.15/6)

print(B)
##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] 0.025 0.025 0.308 0.025 0.025 0.025
## [2,] 0.450 0.025 0.308 0.025 0.025 0.025
## [3,] 0.450 0.025 0.025 0.025 0.025 0.025
## [4,] 0.025 0.025 0.025 0.025 0.450 0.875
## [5,] 0.025 0.025 0.308 0.450 0.025 0.025
## [6,] 0.025 0.025 0.025 0.450 0.450 0.025

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

# uniform rank vector
r <- rep(1/6,6)

# power iterations performed on B until convergence
PwrItr <- function(Bmat, npower, rvec){
  Reduce('%*%', replicate(npower, Bmat, simplify = FALSE)) %*% rvec
}

PwrItr(B,10,r)
##       [,1]
## [1,] 0.018
## [2,] 0.026
## [3,] 0.020
## [4,] 0.169
## [5,] 0.093
## [6,] 0.129
PwrItr(B,20,r)
##       [,1]
## [1,] 0.011
## [2,] 0.016
## [3,] 0.012
## [4,] 0.103
## [5,] 0.057
## [6,] 0.079
PwrItr(B,30,r)
##        [,1]
## [1,] 0.0066
## [2,] 0.0096
## [3,] 0.0074
## [4,] 0.0628
## [5,] 0.0347
## [6,] 0.0480
PwrItr(B,40,r)
##        [,1]
## [1,] 0.0041
## [2,] 0.0059
## [3,] 0.0045
## [4,] 0.0382
## [5,] 0.0211
## [6,] 0.0292

The matrix A set up is not converging. We have noticed that node 2 has no outlinks and the whole row is represented by 0’s. This gives the node 2 zero rank. Let us change each element in this row to 1/6.

A <- matrix(c(0, 1/2, 1/2, 0, 0, 0,
              1/6, 1/6, 1/6, 1/6, 1/6, 1/6,
              1/3, 1/3, 0, 0, 1/3, 0,
              0, 0, 0, 0, 1/2, 1/2,
              0,0, 0, 1/2, 0, 1/2,
              0, 0, 0, 1, 0, 0), nrow=6)


print(A)
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  0.0 0.17 0.33  0.0  0.0    0
## [2,]  0.5 0.17 0.33  0.0  0.0    0
## [3,]  0.5 0.17 0.00  0.0  0.0    0
## [4,]  0.0 0.17 0.00  0.0  0.5    1
## [5,]  0.0 0.17 0.33  0.5  0.0    0
## [6,]  0.0 0.17 0.00  0.5  0.5    0
# Decay value applied and matrix B is formed
d <- 0.85
B <- d * A + (0.15/6)

print(B)
##       [,1] [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] 0.025 0.17 0.308 0.025 0.025 0.025
## [2,] 0.450 0.17 0.308 0.025 0.025 0.025
## [3,] 0.450 0.17 0.025 0.025 0.025 0.025
## [4,] 0.025 0.17 0.025 0.025 0.450 0.875
## [5,] 0.025 0.17 0.308 0.450 0.025 0.025
## [6,] 0.025 0.17 0.025 0.450 0.450 0.025
PwrItr(B,10,r)
##       [,1]
## [1,] 0.052
## [2,] 0.074
## [3,] 0.058
## [4,] 0.348
## [5,] 0.200
## [6,] 0.268
PwrItr(B,20,r)
##       [,1]
## [1,] 0.052
## [2,] 0.074
## [3,] 0.057
## [4,] 0.349
## [5,] 0.200
## [6,] 0.269
PwrItr(B,30,r)
##       [,1]
## [1,] 0.052
## [2,] 0.074
## [3,] 0.057
## [4,] 0.349
## [5,] 0.200
## [6,] 0.269
PwrItr(B,40,r)
##       [,1]
## [1,] 0.052
## [2,] 0.074
## [3,] 0.057
## [4,] 0.349
## [5,] 0.200
## [6,] 0.269

(3) Compute the eigen-decomposition of B and verify that you get an eigenvalue of 1 as 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.

# Get eigenvalues and eigenvectors for matrix B
eigenvalues <- eigen(B)$values
eigenvectors <- eigen(B)$vectors

# All eigenvalues
eigenvalues
## [1]  1.000+0i  0.576+0i -0.425+0i -0.425-0i -0.350+0i -0.085+0i
# Largest eigenvalue
max(as.numeric(eigenvalues))
## Warning: imaginary parts discarded in coercion
## [1] 1

The largest eigenvalue is 1.

# Get the corresponding eigenvector of eigenvalue 1
as.data.frame(eigenvectors[,1])
##   eigenvectors[, 1]
## 1           0.10+0i
## 2           0.15+0i
## 3           0.12+0i
## 4           0.70+0i
## 5           0.40+0i
## 6           0.54+0i
ev <- eigenvectors[,1]

# scale eigenvectors
scale_ev <- data.frame(as.numeric(ev*(1/sum(ev))))
scale_ev
##   as.numeric.ev....1.sum.ev...
## 1                        0.052
## 2                        0.074
## 3                        0.057
## 4                        0.349
## 5                        0.200
## 6                        0.269
# Sum of scaled eigenvectors should be 1
sum(scale_ev)
## [1] 1

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

library(igraph)
## Warning: package 'igraph' was built under R version 3.2.5
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
# plot network
g1 <- graph(c(1,2, 1,3, 3,1, 3,2, 3,5, 4,5, 4,6, 5,6, 5,4, 6,4))
g1$weight <- c(1/2,1/2,1/3,1/3,1/3,1/2,1/2,1/2,1/2,1)

plot(g1)

page.rank(g1)$vector
## [1] 0.052 0.074 0.057 0.349 0.200 0.269

The page.rank vector and the A converged vector are the same!