hwk 10 3 b

gistfile1.txt

power2 <- function(B1 = B1, B2 = B2, v1.old = u1, v2.old = u2, p = .15) {
  N <- length(v1.old)
  B1.p <- matrix(0, ncol=5242, nrow=5242)
  B1.p <- sapply(1:5242, function(x) (1-p)/sum(B1[x,]) * B1[x,])
  B2.p <- (p/N) * B2
  
  v1.new <- B1.p %*% v1.old + B2.p %*% v1.old
  v1.new <- v1.new / norm(v1.new)
  
  v2.new <- B1.p %*% v2.old + B2.p %*% v2.old
  v2.new <- v2.new - as.numeric( (t(v1.new) %*% v2.new) / (t(v1.new) %*% v1.new) ) * v1.new
  v2.new <- v2.new/norm(v2.new)
  
  i <- 1
  
  while( i < 100) {
    v1.old <- v1.new
    v2.old <- v2.new
    
    v1.new <- B1.p %*% v1.old + B2.p %*% v1.old
    v1.new <- v1.new / norm(v1.new)
  
    v2.new <- B1.p %*% v2.old + B2.p %*% v2.old
    v2.new <- v2.new <- v2.new - as.numeric( (t(v1.new) %*% v2.new) / (t(v1.new) %*% v1.new) ) * v1.new
    v2.new <- v2.new/norm(v2.new)
    i <- i + 1
  }
  return(list(vector1 = v1.new, 
              value1 = t(v1.new) %*% (B1.p + B2.p) %*% v1.new, 
              topnode1 = which.max(v1.new),
              vector2 = v2.new, 
              value2 = t(v2.new) %*% (B1.p + B2.p) %*% v2.new, 
              topnode2 = which.max(v2.new),
              iter = i
              ))
}

##results
#> out$value1
#     [,1]
#[1,]    1
#> out$topnode1
#[1] 2775
#> out$value2
#       [,1]
#[1,] 0.2558
$> out$topnode2
#[1] 4784