results from foreach loop in R

Question

I have a function that I need to run on 2000 data frames. Each iteration is taking a very long time i.e almost 40 minutes and hence I'm using the 'foreach' package in R. I have generated the data in the following way:

library(foreach)
library(doParallel)
library(data.table)
library(matrixStats)
# DATA
datalist <- list()
for (i in 1:2000){
  set.seed(i)
  x1 <- rnorm(600,0.05,0.3)
  x2 <- rnorm(600,-2,0.25)
  data_2 <- data.frame(x1,x2)
  lin_pred <- 1+(0.2 * data_2[1]) + (1.2*data_2[2]) 
  prob <- 1/(1+exp(-lin_pred))
  index <- rep(1:100, each = 6)
  data <- data.frame(index,prob)
  data_split <- split(data,f=data$index)
  create_y <- function(p){
    
    y <- c()
    y_1 <- rbinom(1,1,p[1,2])
    y <- append(y,y_1)
    for(i in 2:6){
      pr <- p[i,2] + 0.05*(y[i-1]-p[i-1,2])
      u <- rbinom(1,1,pr)
      y <- append(y,u)
    }
    return(y)
  }
  
  res <- lapply(data_split,create_y)
  y <- data.frame(x=unlist(res))
  final_data <- data.frame(index,y,x1,x2)
  datalist[[i]] <- final_data
}

Now I have defined my objective function which will be optimized for each of the data set in the list defined above:

## Objective function:
dpd_tdependent <- function(x, fixed = c(rep(FALSE,5))){
  params <- fixed
  dpd <- function(p){
    params[!fixed] <- p
    alpha <- params[1]
    beta_0 <- params[2]
    beta_1 <- params[3]
    beta_2 <- params[4]
    rho <- params[5]
    add_pi <- function(d){
      k <- beta_0+(d[3]*beta_1)+(d[4]*beta_2)
      k1 <- exp(k)/(1+exp(k))
      p <- c()
      p <- append(p,k1[1,1])
      for(i in 2:6){
        u <- k1[i,1] + rho*(d[i-1,2]-k1[i-1,1])
        p <- append(p,u)
      }
      d <- cbind(d,p)
    }
    dat_split <- split(x , f  = x$index)
    result <- lapply(dat_split, add_pi)
    
    result <- rbindlist(result)
    result <- as.data.frame(result)
    colnames(result) <- c('index','y','x1','x2','exp_prob')
    result_split <- split(result, f = result$index)
    ## First expression
    full_prob <- function(d){
      k <- as.data.frame(expand.grid(c(0,1),c(0,1),c(0,1),c(0,1),c(0,1),c(0,1)))
      k <- as.data.frame(t(k))
      val <- c()
      for(j in 1:ncol(k)){
        d[2]<- k[j]
        m <- d[1,5]^d[1,2] * ((1-d[1,5])^(1-d[1,2]))
        for(i in 2:nrow(d)){
          c <- 1+ (rho*((d[i,2]-d[i,5])*(d[i-1,2]-d[i-1,5]))/(sqrt(abs(d[i,5]*d[i-1,5]*(1-d[i,5])*(1-d[i-1,5])))))
          m <- m* d[i,5]^d[i,2] * (1-d[i,5])^(1-d[i,2]) *c 
        }
        val <- append(val,m)
      }
      val <- val^(1+alpha)
      return(sum(val))
    }
    first_exp <- lapply(result_split,full_prob)
    first_exp <- as.vector(unlist(first_exp))
    
    ## Second expression:
    compute_prob <- function(d){
      m <- d[1,5]^d[1,2] * ((1-d[1,5])^(1-d[1,2]))
      for(i in 2:nrow(d)){
        c <- 1+ (rho*((d[i,2]-d[i,5])*(d[i-1,2]-d[i-1,5]))/(sqrt(abs(d[i,5]*d[i-1,5]*(1-d[i,5])*(1-d[i-1,5])))))
        m <- m* d[i,5]^d[i,2] * (1-d[i,5])^(1-d[i,2]) *c 
      }
      
      return(m^alpha)
    }
    second_exp <- lapply(result_split,compute_prob)
    second_exp <- as.vector(unlist(second_exp))
    
    final_res <- first_exp - ((1+1/alpha)*(second_exp))
    final_result <- sum(final_res)
    
  }
}

Lastly, I'm using the foreach package to hasten the process:

cl = makeCluster(6)
registerDoParallel(cl)




mse <- matrix(,nrow=2000,ncol=5)
foreach(i = 1:2000, .packages = c('data.table','matrixStats'), .export = c('mse','datalist'))%dopar%{
  beta <- rbind(1,0.2,1.2,0.05)
  val <- dpd_tdependent(datalist[[i]], c(0.7,FALSE,FALSE,FALSE,FALSE))
  b_s <- as.vector(optim(c(beta_0 =0.7, beta_1 =0.05 ,beta_2 = 0.9,rho=0.001),val)$par)
  conv <- optim(c(beta_0 =0.7, beta_1 =0.05 ,beta_2 = 0.9,rho = 0.001),val)$convergence
  k <- (b_s -beta)
  k <- append(k,conv)
  mse[i,] <- rbind(k)
  if(colCounts(mse , value = 0, na.rm = TRUE)[5] == 500) break
}
stopCluster(cl)

Now, my problem is that the matrix mse is not being populated with the desired values after running the optim function inside the foreach package. I am aware that I can use the .combine feature inside the foreach function to return a matrix, but how will I assign a name to such a matrix? I mean if I cannot assign a name then, how will I check the last condition(inside if) and break?

I truly appreciate any help in this regard. Thank you.

Answer 1

I don't believe you should be trying to modify a global variable from within each worker. See my comment above and link. You shouldn't be checking within iteration process if 500 iterations have convergence=0, because that information is not available to each iteration. The below is one option to return what you want

cl = makeCluster(6)
registerDoParallel(cl)

mse = foreach(i = 1:2000, .packages = c('data.table','matrixStats')) %dopar%{
  beta <- rbind(1,0.2,1.2,0.05)
  val <- dpd_tdependent(datalist[[i]], c(0.7,FALSE,FALSE,FALSE,FALSE))
  optim_sol <- optim(c(beta_0 =0.7, beta_1 =0.05 ,beta_2 = 0.9,rho=0.001),val)
  b_s <- optim_sol$par
  conv <- optim_sol$convergence
  c(b_s-beta,conv,i)
}
mse <- matrix(unlist(m),nrow=2000, byrow=T)

stopCluster(cl)