hwk10

gistfile1.txt

Elaine Ayo

Math 504 / Homework 10

March 29, 2015

```{r echo=FALSE, message=FALSE}
library(MASS)
library(Matrix)
library(ggplot2)
library(reshape2) 
library(pryr)
library(microbenchmark)
library(parallel)
```


2. The resulting function calculates a dominant eigenvector/value roughly equal to R's eigen function, but the signs are flipped on the vectors. 
```{r}
norm <- function(x) sqrt(sum(x^2))

set.seed(693)
M <- matrix(runif(9), ncol=3, nrow=3)
v <- rep(1,3)

eig.pm <- function(M = M, v.old = v, tol = 10^-10) {
  v.new <- M %*% v.old
  v.new <- v.new / norm(v.new)
  i <- 0
  
  while( norm(v.new - v.old) > tol) {
    v.old <- v.new
    v.new <- M %*% v.old
    v.new <- v.new / norm(v.new)
    i <- i + 1
  }
  return(list(pm.vector = v.new, 
              pm.value = t(v.new)%*%M%*%v.new, 
              R.vector = eigen(M)$vectors[,1], 
              R.value = max(eigen(M)$values),
              matrix = M
              ))
}

eig.pm(M = M, v.old = v, tol = 10^-10)
```


3. a.

```{r}
#import data
data <- read.delim("~/Dropbox/numericalmethods/ca-GrQc.txt", header=TRUE)

same <- which(data$SNode == data$ENode)

#there are 12 nodes that link to themselves

#create sparse matrix B1
#find all unique nodes, there are 5,242 unique numbers as expected
nodes <- unique(c(data$SNode, data$ENode))
nodes <- nodes[order(nodes)]
trans <- data.frame(id=1:length(nodes), num = nodes)

#transform the numbering in the data to reflect sequential numbers to make generating the matrix easier
data2<- merge(data, trans, by.x="SNode", by.y="num", all.x=TRUE)
names(data2)[3]<-"i"
data2<-merge(data2,trans,by.x="ENode",by.y="num",all.x=TRUE)
names(data2)[4]<-"j"

data2 <- data2[,3:4]

#construct B1
B1 <- matrix(0, nrow=5242, ncol=5242)

for (x in 1:28980) {
  B1[data2$i[x], data2$j[x]] <- 1
}

#check all edges accounted for, there should be 28,980
sum(B1)

#construct B2
B2 <- matrix(1, nrow=5242, ncol=5242)

#implement modified power method
norm <- function(x) sqrt(sum(x^2))

set.seed(693)
v <- runif(5242)

power <- function(B1 = B1, B2 = B2, v.old = v, tol = 10^-10, p = .15) {
  N <- length(v)
  v.new <- (1-p) * B1 %*% v.old + (p/N) * B2 %*% v.old
  v.new <- v.new / norm(v.new)
  i <- 0
  
  while( norm(v.new - v.old) > tol) {
    v.old <- v.new
    v.new <- (1-p) * B1 %*% v.old + (p/N) * B2 %*% v.old
    v.new <- v.new / norm(v.new)
    i <- i + 1
  }
  return(list(vector = v.new, 
              value = t(v.new) %*% ((1-p) * B1) %*% v.new + t(v.new) %*% ((p/N) * B2) %*% v.new, 
              topnode = which.max(v.new),
              iter = i
              ))
}

full.time <- list()
results <- list()

#now test a lot of p's 
p <- seq(0, 1, .15)
p[8] <- 1


full.time <- lapply(1:2, function(x) system.time(results[[x]] <- power(B1 = B1, B2 = B2, v.old = v, tol = 10^-10, p = x)) )

#use parallel processing
full.time <- mclapply(2:3, function(x) system.time(results[[x]] <- power(B1 = B1, B2 = B2, v.old = v, tol = 10^-10, p = x)), mc.silent = FALSE, mc.cores = getOption("cores",4), mc.cleanup = TRUE)

#now run power using sparse matrices
B1.sparse <- Matrix(0, nrow=5242, ncol=5242, sparse=TRUE)

for (x in 1:28980) {
  B1.sparse[data2$i[x], data2$j[x]] <- 1
}

#check all edges accounted for, there should be 28,980
sum(B1.sparse)

full.time.sp <- list()
results.sp <- list()

full.time.sp <- mclapply(2:3, function(x) system.time(results.sp[[x]] <- power(B1 = B1.sparse, B2 = B2, v.old = v, tol = 10^-10, p = x)), mc.silent = FALSE, mc.cores = getOption("cores",4), mc.cleanup = TRUE)
```