Data 605 HW10

•Form the A matrix. Then, introduce decay and form the B matrix

# Transition matrix A
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)

# Take Decay d = 0.85

d <- 0.85

B <- 0.85*A + (0.15/6)

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

r <- matrix(c(1/6,1/6,1/6,1/6,1/6,1/6),nrow=6)

# perform power iterations on B
iterations <- function(M, r, n) {
  Bn = diag(nrow(M)) 
  # Caculate B n times
  for (i in 1:n)
  {
    Bn = Bn %*% M
  }
  return (Bn %*% r)
}

# Different values of n, to see convergence
eig_vec1<-iterations(B, r, 30) #Converges
#iterations(B, r, 40)
eig_vec1

##            [,1]
## [1,] 0.05170475
## [2,] 0.07367927
## [3,] 0.05741242
## [4,] 0.34870367
## [5,] 0.19990381
## [6,] 0.26859607

• Compute the eigen-decomposition of B and verify that you indeed get an eigenvalueof 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.

# Decomposing B
eigen_decomp <- eigen(B)

# Maximum eigen value we get is indeed 1
max_value<-which.max(eigen_decomp$values)

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

# check this eigenvector has all positive entries and it sums to 1
eig_vec2 <- as.numeric((1/sum(eigen_decomp$vectors[,1]))*eigen_decomp$vectors[,1])
sum(eig_vec2)

## [1] 1

• 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)

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

graph_A <- graph.adjacency(t(A), weighted = TRUE, mode = "directed")
plot(graph_A)

eig_vec3 <- page.rank(graph_A)$vector
eig_vec3

## [1] 0.05170475 0.07367926 0.05741241 0.34870369 0.19990381 0.26859608

Verify that you do get the same PageRank vector as the two approaches above.

results <- cbind(eig_vec1, eig_vec2, eig_vec3)
colnames(results) <- c('method1', 'method2', 'method3')
results

##         method1    method2    method3
## [1,] 0.05170475 0.05170475 0.05170475
## [2,] 0.07367927 0.07367926 0.07367926
## [3,] 0.05741242 0.05741241 0.05741241
## [4,] 0.34870367 0.34870369 0.34870369
## [5,] 0.19990381 0.19990381 0.19990381
## [6,] 0.26859607 0.26859608 0.26859608