## The matrix is not in its proper form we have to transpose it so that the row elements correspond to its inlink and column wise it corresponds to its outlink each column represents a page. ## I also introduced 1/6 in the second row so that the matrix behaves well..
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)
## create the decay so

d <- 0.85

B <- (d* A) + (0.15/6)
## create the r vector
r <- matrix(c(1/6,1/6,1/6,1/6,1/6,1/6),nrow=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

pow_iter <- function(A,r,n){
  B = diag(nrow(A))
  for (i in 1:n){
    ## keep updating the B matrix by number of iterations
    B = B %*% A
  
  }
  return (B %*% r)
}  
## experiment with numbers of iterations.. 
pow_iter(B,r,1)
##            [,1]
## [1,] 0.09583333
## [2,] 0.16666667
## [3,] 0.11944444
## [4,] 0.26111111
## [5,] 0.16666667
## [6,] 0.19027778
pow_iter(B,r,10)
##            [,1]
## [1,] 0.05205661
## [2,] 0.07428990
## [3,] 0.05782138
## [4,] 0.34797267
## [5,] 0.19975859
## [6,] 0.26810085
pow_iter(B,r,14)
##            [,1]
## [1,] 0.05174349
## [2,] 0.07374658
## [3,] 0.05745753
## [4,] 0.34862496
## [5,] 0.19988600
## [6,] 0.26854144
pow_iter(B,r,20)
##            [,1]
## [1,] 0.05170616
## [2,] 0.07368173
## [3,] 0.05741406
## [4,] 0.34870083
## [5,] 0.19990313
## [6,] 0.26859408
(vecto1 <- pow_iter(B,r,30))
##            [,1]
## [1,] 0.05170475
## [2,] 0.07367927
## [3,] 0.05741242
## [4,] 0.34870367
## [5,] 0.19990381
## [6,] 0.26859607
pow_iter(B,r,40)
##            [,1]
## [1,] 0.05170475
## [2,] 0.07367926
## [3,] 0.05741241
## [4,] 0.34870369
## [5,] 0.19990381
## [6,] 0.26859608

It seems that the B vector converges after 15 to 30 iterations since we are now getting approximately the same values:

Compute the eigen-decomposition of B and verify that you indeed get an eigenvalue of 1

eigen(B)
## 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
## verify that we got the max eigenvalue of 1:
eigenvalue = eigen(B)$values
max(Re(eigenvalue[abs(Im(eigenvalue)) < 1e-6]))
## [1] 1
## We can see that the first set of eigenvectors are all positive as given by the theorem where all the vectors are positive..
eigen(B)$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
## The eigen function is normalizing the eigenvector so we have to find the cancel the normalization by  
vecto2 <- as.numeric((1/sum(eigen(B)$vectors[,1]))*eigen(B)$vectors[,1])
vecto2
## [1] 0.05170475 0.07367926 0.05741241 0.34870369 0.19990381 0.26859608

###Using the igraph function:

library(igraph)
## Warning: package 'igraph' was built under R version 4.1.3
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
## I didn't how to use this... nned to use the transpose for the argument..and plot it to use for page rank method
graphA <-graph.adjacency(t(A),weighted = TRUE,mode="directed")
plot(graphA)

vecto3 <- page.rank(graphA)$vector
vecto3
## [1] 0.05170475 0.07367926 0.05741241 0.34870369 0.19990381 0.26859608
data <- cbind(vecto1,vecto2,vecto3)
colnames(data) <- c("Iteration Method","Eigen Function","Page Rank Function")
data
##      Iteration Method Eigen Function Page Rank Function
## [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