page rank solution -problem 1

library(tidyverse)

## Warning: package 'tibble' was built under R version 4.0.4

library(openintro)
library(Matrix)
if (!require('matrixcalc')) install.packages('matrixcalc')

## Warning: package 'matrixcalc' was built under R version 4.0.5

if (!require('igraph')) install.packages('igraph')

## Warning: package 'igraph' was built under R version 4.0.4

Exercise 1

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

p1 <- c(0, 1/2, 1/2, 0, 0, 0)
p2 <- rep(1/6, 6) # we adjust this from 0 so that we have an equal prob to land on any page;   if we leave at 0 (indicating there are no outgoing links from p2), our page rank will not converge because of a 'dangling node'
p3 <- c(1/3, 1/3, 0, 0, 1/3, 0)
p4 <- c(0, 0, 0, 0, 1/2, 1/2)
p5 <- c(0, 0, 0, 1/2, 0, 1/2)
p6 <- c(0, 0, 0, 1, 0, 0)
A <- matrix(c(p1, p2, p3, p4, p5, p6), 6)

#confirm that the total probability for each col = 1
colSums(A)

## [1] 1 1 1 1 1 1

#Introduce decay--.85 is the damping factor
B <- 0.85 * A + 0.15/nrow(A)

Exercise 2

Start with a uniform rank vector r and perform power iterations on B till convergence. That is, compute the solution r = B^n × r. Attempt this for a sufficiently large n so that r actually converges. (5 Points)

The following function is created to perform power iterations on B until convergence, utilizing a uniform rank vector \[r^T = \left[ \begin{array}{c} \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} \end{array} \right]\]

r <- rep(1/nrow(A), nrow(A))

b<-cbind(matrix.power(B, 10) %*% r,
      matrix.power(B, 20) %*% r,
      matrix.power(B, 30) %*% r,
      matrix.power(B, 40) %*% r,
      matrix.power(B, 50) %*% r,
      matrix.power(B, 60) %*% r)

b

##            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## [1,] 0.05205661 0.05170616 0.05170475 0.05170475 0.05170475 0.05170475
## [2,] 0.07428990 0.07368173 0.07367927 0.07367926 0.07367926 0.07367926
## [3,] 0.05782138 0.05741406 0.05741242 0.05741241 0.05741241 0.05741241
## [4,] 0.34797267 0.34870083 0.34870367 0.34870369 0.34870369 0.34870369
## [5,] 0.19975859 0.19990313 0.19990381 0.19990381 0.19990381 0.19990381
## [6,] 0.26810085 0.26859408 0.26859607 0.26859608 0.26859608 0.26859608

PR1 <- matrix.power(B, 40) %*% r # it looks like convergence  occurs at 40 iterations
PR1

##            [,1]
## [1,] 0.05170475
## [2,] 0.07367926
## [3,] 0.05741241
## [4,] 0.34870369
## [5,] 0.19990381
## [6,] 0.26859608

Exercise 3

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.(10 points)

decomp <- eigen(B)
decomp

## eigen() decomposition
## $values
## [1]  1.00000000+0i  0.57619235+0i -0.42500000+0i -0.42500000-0i -0.34991524+0i
## [6] -0.08461044+0i
## 
## $vectors
##              [,1]          [,2]                        [,3]
## [1,] 0.1044385+0i  0.2931457+0i  2.945054e-15+5.507002e-22i
## [2,] 0.1488249+0i  0.5093703+0i -1.223015e-15-0.000000e+00i
## [3,] 0.1159674+0i  0.3414619+0i -2.241513e-15-6.032865e-22i
## [4,] 0.7043472+0i -0.5890805+0i -7.071068e-01+0.000000e+00i
## [5,] 0.4037861+0i -0.1413606+0i  7.071068e-01+0.000000e+00i
## [6,] 0.5425377+0i -0.4135367+0i  0.000000e+00-2.145851e-08i
##                             [,4]           [,5]            [,6]
## [1,]  2.945054e-15-5.507002e-22i -0.06471710+0i -0.212296003+0i
## [2,] -1.223015e-15+0.000000e+00i  0.01388698+0i  0.854071294+0i
## [3,] -2.241513e-15+6.032865e-22i  0.07298180+0i -0.363638739+0i
## [4,] -7.071068e-01+0.000000e+00i -0.66058664+0i  0.018399984+0i
## [5,]  7.071068e-01-0.000000e+00i  0.73761812+0i -0.304719509+0i
## [6,]  0.000000e+00+2.145851e-08i -0.09918316+0i  0.008182973+0i

Re(eigen(B)$values[1])

## [1] 1

# Sum of the vector is greater than 1 for the eigen value 1
#sum(eigen(B)$vectors[,1])

PR2 <- as.numeric(decomp$vectors[,which.max(decomp$values)]) #get vectors associated with largest eigenvalue == 1

## Warning in which.max(decomp$values): imaginary parts discarded in coercion

sum(PR2)

## [1] 2.019902

# Change it to unit vector
PR2 <- (1/sum(PR2))*PR2 #normalize
PR2

## [1] 0.05170475 0.07367926 0.05741241 0.34870369 0.19990381 0.26859608

sum(PR2)

## [1] 1

# Difference between power method and eigen vector is negligible
sum(PR1 - PR2)

## [1] 3.469447e-17

Exercise 4

Use the graph package in R and its page.rank method to compute the Page Rank of the graph as given in A. Note that you don’t need to apply decay. The package starts with a connected graph and applies decay internally. Verify that you do get the same PageRank vector as the two approaches above(10 points).

#Converting to graph from adjacency matrix
g1 = graph_from_adjacency_matrix(t(A),weighted = T)

#Plot the graph
plot(g1)

#Resultant vector
page_rank(g1)$vector

## [1] 0.05170475 0.07367926 0.05741241 0.34870369 0.19990381 0.26859608

#we can conclude that we are getting the same results for eigen,graph and power iteration