Computational Mathematics - PageRank

##Introduction You’ll verify for yourself that PageRank works by performing calculations on a small universe of web pages.

Let’s use the 6 page universe that we had in the previous discussion For this directed graph, perform the following calculations in R.

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

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.

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)

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)

library(kableExtra)
library(matrixcalc)
library(igraph)

## 
## Attaching package: 'igraph'

## The following object is masked from 'package:matrixcalc':
## 
##     %s%

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

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

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,1/2,1/2,0,0,0,1/2,0,1/2,0,0,0,1,0,0), nrow = 6, byrow = TRUE)

#kable(A, digits = 2, caption="Matrix A")

n <- 6
B <- 0.85*A + 0.15/n

kable(B, digits = 2, caption = "Matrix B")

r <- matrix(c(.167,.167,.167,.167,.167,.167))
kable(r, caption = "Rank Vector") 

r2 <- B %*% r
r3 <- B %*% r2
r3

##            [,1]
## [1,] 0.10312250
## [2,] 0.02150125
## [3,] 0.12323208
## [4,] 0.16345125
## [5,] 0.16345125
## [6,] 0.16345125

r2 <- A %*% r
r2

##       [,1]
## [1,] 0.167
## [2,] 0.000
## [3,] 0.167
## [4,] 0.167
## [5,] 0.167
## [6,] 0.167

r3 <- A %*% r2
r3

##           [,1]
## [1,] 0.0835000
## [2,] 0.0000000
## [3,] 0.1113333
## [4,] 0.1670000
## [5,] 0.1670000
## [6,] 0.1670000

r4 <- A %*% r3
r4

##            [,1]
## [1,] 0.05566667
## [2,] 0.00000000
## [3,] 0.08350000
## [4,] 0.16700000
## [5,] 0.16700000
## [6,] 0.16700000

r5 <- A %*% r4

matrix.power(A,2)%*%r

##           [,1]
## [1,] 0.0835000
## [2,] 0.0000000
## [3,] 0.1113333
## [4,] 0.1670000
## [5,] 0.1670000
## [6,] 0.1670000

matrix.power(A,3)%*%r

##            [,1]
## [1,] 0.05566667
## [2,] 0.00000000
## [3,] 0.08350000
## [4,] 0.16700000
## [5,] 0.16700000
## [6,] 0.16700000

matrix.power(A,4)%*%r

##            [,1]
## [1,] 0.04175000
## [2,] 0.00000000
## [3,] 0.07422222
## [4,] 0.16700000
## [5,] 0.16700000
## [6,] 0.16700000

matrix.power(A,10)%*%r

##            [,1]
## [1,] 0.03343866
## [2,] 0.00000000
## [3,] 0.06683436
## [4,] 0.16700000
## [5,] 0.16700000
## [6,] 0.16700000

matrix.power(A,20)%*%r

##        [,1]
## [1,] 0.0334
## [2,] 0.0000
## [3,] 0.0668
## [4,] 0.1670
## [5,] 0.1670
## [6,] 0.1670

matrix.power(B,2)%*%r

##            [,1]
## [1,] 0.10312250
## [2,] 0.02150125
## [3,] 0.12323208
## [4,] 0.16345125
## [5,] 0.16345125
## [6,] 0.16345125

matrix.power(B,3)%*%r

##            [,1]
## [1,] 0.07996691
## [2,] 0.01845524
## [3,] 0.10007649
## [4,] 0.15738880
## [5,] 0.15738880
## [6,] 0.15738880

matrix.power(B,4)%*%r

##            [,1]
## [1,] 0.06714261
## [2,] 0.01676663
## [3,] 0.08924639
## [4,] 0.15054711
## [5,] 0.15054711
## [6,] 0.15054711

matrix.power(B,10)%*%r

##            [,1]
## [1,] 0.04527118
## [2,] 0.01199487
## [3,] 0.06248074
## [4,] 0.11226120
## [5,] 0.11226120
## [6,] 0.11226120

matrix.power(B,20)%*%r

##             [,1]
## [1,] 0.027561208
## [2,] 0.007305781
## [3,] 0.038049817
## [4,] 0.068395282
## [5,] 0.068395282
## [6,] 0.068395282

matrix.power(B,30)%*%r

##             [,1]
## [1,] 0.016791252
## [2,] 0.004450938
## [3,] 0.023181282
## [4,] 0.041668806
## [5,] 0.041668806
## [6,] 0.041668806

matrix.power(B,40)%*%r

##             [,1]
## [1,] 0.010229820
## [2,] 0.002711668
## [3,] 0.014122850
## [4,] 0.025386099
## [5,] 0.025386099
## [6,] 0.025386099

matrix.power(B,80)%*%r

##              [,1]
## [1,] 0.0014093135
## [2,] 0.0003735735
## [3,] 0.0019456377
## [4,] 0.0034973218
## [5,] 0.0034973218
## [6,] 0.0034973218

matrix.power(B,100)%*%r

##              [,1]
## [1,] 0.0005230912
## [2,] 0.0001386583
## [3,] 0.0007221573
## [4,] 0.0012980918
## [5,] 0.0012980918
## [6,] 0.0012980918

matrix.power(B,200)%*%r

##              [,1]
## [1,] 3.684907e-06
## [2,] 9.767760e-07
## [3,] 5.087224e-06
## [4,] 9.144384e-06
## [5,] 9.144384e-06
## [6,] 9.144384e-06

matrix.power(B,400)%*%r

##              [,1]
## [1,] 1.828625e-10
## [2,] 4.847224e-11
## [3,] 2.524521e-10
## [4,] 4.537876e-10
## [5,] 4.537876e-10
## [6,] 4.537876e-10

d <- 0.85
n <- 6
(matrix.power(A,2)%*%r)*d + (1-d)/n

##           [,1]
## [1,] 0.0959750
## [2,] 0.0250000
## [3,] 0.1196333
## [4,] 0.1669500
## [5,] 0.1669500
## [6,] 0.1669500

(matrix.power(A,4)%*%r)*d + (1-d)/n

##            [,1]
## [1,] 0.06048750
## [2,] 0.02500000
## [3,] 0.08808889
## [4,] 0.16695000
## [5,] 0.16695000
## [6,] 0.16695000

(matrix.power(A,10)%*%r)*d + (1-d)/n

##            [,1]
## [1,] 0.05342286
## [2,] 0.02500000
## [3,] 0.08180921
## [4,] 0.16695000
## [5,] 0.16695000
## [6,] 0.16695000

ev <- eigen(A)
ev$values

## [1]  1.0000000 -0.5000000 -0.5000000  0.4082483 -0.4082483  0.0000000

ev$vectors

##           [,1]       [,2]       [,3]      [,4]       [,5]       [,6]
## [1,] 0.1118034 -0.5345225 -0.5345225 0.7745967  0.7745967 -0.5773503
## [2,] 0.0000000  0.0000000  0.0000000 0.0000000  0.0000000  0.5773503
## [3,] 0.2236068  0.5345225  0.5345225 0.6324555 -0.6324555 -0.5773503
## [4,] 0.5590170 -0.2672612 -0.2672612 0.0000000  0.0000000  0.0000000
## [5,] 0.5590170 -0.2672612 -0.2672612 0.0000000  0.0000000  0.0000000
## [6,] 0.5590170  0.5345225  0.5345225 0.0000000  0.0000000  0.0000000

#Error in this matrix further investigation required. 
#matrix.inverse(ev$vectors)

evB <- eigen(B)
evB$values

## [1]  0.95165271+0i -0.42500000+0i -0.42500000-0i  0.41436582+0i -0.34733611+0i
## [6] -0.01868242+0i

evB$vectors

##               [,1]                        [,2]                        [,3]
## [1,] -0.2159123+0i -5.345225e-01-0.000000e+00i -5.345225e-01+0.000000e+00i
## [2,] -0.0572329+0i -5.200810e-17+3.664536e-10i -5.200810e-17-3.664536e-10i
## [3,] -0.2980792+0i  5.345225e-01+0.000000e+00i  5.345225e-01+0.000000e+00i
## [4,] -0.5358032+0i -2.672612e-01+0.000000e+00i -2.672612e-01-0.000000e+00i
## [5,] -0.5358032+0i -2.672612e-01+0.000000e+00i -2.672612e-01-0.000000e+00i
## [6,] -0.5358032+0i  5.345225e-01-0.000000e+00i  5.345225e-01+0.000000e+00i
##                [,4]            [,5]           [,6]
## [1,] -0.77869913+0i -0.775079691+0i  0.55060201+0i
## [2,] -0.07537631+0i  0.009090375+0i -0.61511270+0i
## [3,] -0.61034825+0i  0.631781588+0i  0.56386946+0i
## [4,]  0.07169632+0i  0.002637034+0i -0.01322899+0i
## [5,]  0.07169632+0i  0.002637034+0i -0.01322899+0i
## [6,]  0.07169632+0i  0.002637034+0i -0.01322899+0i

matrix.inverse(evB$vectors)

##                             [,1]                        [,2]
## [1,] -8.097330e-02+0.000000e+00i -1.171353e-01+0.000000e+00i
## [2,]  0.000000e+00-3.368001e-09i  0.000000e+00-2.020801e-08i
## [3,]  0.000000e+00+3.368001e-09i  0.000000e+00+2.020801e-08i
## [4,] -5.688106e-01-0.000000e+00i -1.152219e+00-0.000000e+00i
## [5,] -6.481031e-01+0.000000e+00i  1.449150e-01-0.000000e+00i
## [6,]  6.765863e-02+0.000000e+00i -1.471484e+00+0.000000e+00i
##                             [,3]                        [,4]
## [1,] -9.026184e-02+0.000000e+00i -7.644771e-01+0.000000e+00i
## [2,]  0.000000e+00-6.736003e-09i  0.000000e+00-7.181207e+07i
## [3,]  0.000000e+00+6.736003e-09i  0.000000e+00+7.181207e+07i
## [4,] -6.843066e-01-0.000000e+00i  4.840239e-01-0.000000e+00i
## [5,]  7.864378e-01-0.000000e+00i -6.603705e+00-0.000000e+00i
## [6,]  1.038760e-01+0.000000e+00i -2.101631e-04+2.047285e-10i
##                       [,5]           [,6]
## [1,] -0.42238228+       0i -0.58414116+0i
## [2,]  0.00000000+71812067i  0.62360956+0i
## [3,]  0.00000000-71812067i  0.62360956-0i
## [4,] -0.07236477+       0i  0.32132564+0i
## [5,]  7.43218116+       0i -1.02030281-0i
## [6,]  0.07243956-       0i -0.00010272-0i

PageR_Mat <- 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,1/2,1/2,0,0,0,1/2,0,1/2,0,0,0,1,0,0), nrow = 6, byrow = TRUE)
network <- graph_from_adjacency_matrix(PageR_Mat)
plot(network)