library(igraph)         # for graph_from_data_frame and page_Rank





Form the Matrix

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


First create a function that performs power iteration until a convergence threshold is reached.


pow_it_until<-function(A, r, it) {
  
  epsilon<-.00001
  
  r_t<-r   # intialize
  
  for ( i in 1:it) {
    r_t_plus1<- r_t  %*% t(A)

    if (sum(abs(r_t_plus1-r_t)) < epsilon) {
      print(paste("Breaking after", i, " iterations"))
      break
      
    }
    r_t<-r_t_plus1
  }
  
  return(r_t_plus1)
}


Now create a function that defines a decay matrix.


create_decay_matrix<-function(A, d) {
  
  return(decay * t(A) + (1-decay)/n)
}


This is the matrix…


\[A = \begin{vmatrix} 0&\frac{1}{2}&\frac{1}{2}&0&0&0 \\ 0&0&0&0&0&0\\ \frac{1}{3}&\frac{1}{3}&0&0&\frac{1}{3}&0 \\ 0&0&0&0&\frac{1}{2}&\frac{1}{2}\\ 0&0&0&\frac{1}{2}&0&\frac{1}{2}\\ 0&0&0&1&0&0\\ \end{vmatrix}\]


A <- matrix(c(0,.5,.5,0,0,0,  0,0,0,0,0,0, 1/3, 1/3 ,0,0, 1/3, 0, 0,0,0,0,.5, .5, 0,0,0, .5, 0, .5, 0,0,0, 1,0,0 ), nrow=6, byrow=TRUE)
n<-nrow(A)
urv<-rep(1 / n, n)

# A[2,]<-urv

decay<-.95
max_iterations<-100          # in case it doesnt converge

B<-create_decay_matrix(A,decay)


Power Iteration Vector


Start with a uniform rank vector r and perform power iterations on B till conver- gence. That is, compute the solution r = \(B^n * r\). Attempt this for a sufficiently large n so that r actually converges.

pow_vec<-pow_it_until(B,urv, max_iterations)

print(sprintf("Vector from Iteration :   %s", paste(round(pow_vec,4), sep = " ", collapse = '  ')))
## [1] "Vector from Iteration :   0.0013  0.0019  0.0015  0.0409  0.0211  0.0308"


Eigen Vector


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.


The maximum eigen value will be the first one.


The given vector is just one of infinite vectors. To equate it to above, we need to scale it.


eig<-eigen(B)

eig_vec<-eig$vectors[,1]
eig_val<-eig$values[1]
scalar<-pow_vec[1]/eig_vec[1]
eig_vec<-scalar * eig_vec

print(sprintf("Vector from eigen() :   %s", paste(round(eig_vec,4), sep = " ", collapse = '  ')))

print(sprintf("Value from eigen() :   %.5f",eig_val))
## [1] "Vector from eigen() :   0.0013  0.0019  0.0015  0.0409  0.0211  0.0308"
## [1] "Value from eigen() :   0.98123"


Note: The eigen value isnt 1. I notice with the original matrix its 1. Then as decay increases the eigen value decreases indicating to me a destabilization of the convergence assumption.


In all cases the power vectors and the iteration vectors and the eigen values are correct for the vector…


B_times_eig_vec<-B %*% eig_vec

print(sprintf("B * eig_vec :   %s", paste(round(B_times_eig_vec,4), sep = " ", collapse = '  ')))

print(sprintf("eig_val * eig_vec :   %s", paste(round(eig_val * eig_vec,4), sep = " ", collapse = '  ')))
## [1] "B * eig_vec :   0.0013  0.0019  0.0014  0.0401  0.0207  0.0303"
## [1] "eig_val * eig_vec :   0.0013  0.0019  0.0014  0.0401  0.0207  0.0303"


igraph Vector


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.


links=cbind('from'=c('1', '1', '3', '3','3', '4','4', '5','5','6'),
            'to'=c('2', '3', '1', '2','5', '5','6', '4','6','4')
            )

g<-graph_from_data_frame(links,vertices = c('1','2','3','4','5','6'),directed = T)
plot(g)

g_page_rank<-page_rank(g, directed = T)


igraph_vector<-g_page_rank$vector

scalar<-pow_vec[1]/igraph_vector[1]



print(sprintf("igraph_vec :   %s", paste(round(scalar * igraph_vector,4), sep = " ", collapse = '  ')))
## [1] "igraph_vec :   0.0013  0.0019  0.0014  0.0088  0.005  0.0067"


The first 3 match then the last 3 deviate a little. Im curious why…