CODEWITHSUNDEEP
rmatrixparallel-processing
wolfsatthedoor
wolfsatthedoor

Reputation: 7313

a way to make this array/matrix multiplication faster? in R

I cannot think of a way to make this code faster. Is there an apply function that would run faster? Now, I am using a for each loop to run this loop in parallel, and it still takes a VERY long time. 

ndraws=20000
nhousehold=18831
m=12


    elasticitydraws = array(0,c(m,ndraws,nhousehold))
    MAPelasticity = matrix(0,nhousehold,m)
    medianpricemat = matrix(rnorm(12,15,1),12,1)                                 


# dim(out$betadraw) = 18831, 12, 20000
# dim(medianpricemat) = 12, 1

    library(foreach)
    library(doMC)
    registerDoMC(10)

    elasticitylist = foreach(i=1:nhousehold) %dopar% {

            pricedraws = out$betadraw[i,12,] 
            elasticitydraws[,,i]= probarray[,,i] %*% diag(pricedraws)
            elasticitydraws[,,i] = elasticitydraws[,,i] * as.vector(medianpricemat)
            MAPelasticity[i,] = apply(elasticitydraws[,,i],1,mean)

    }

Upvotes: 1

Views: 1380

Answers (1)

klash
klash

Reputation: 341

The bottleneck in the code is the creation a large dense diagonal matrix and matrix multiplication with this. It is better using a sparse matrix and the Matrix package. This saves memory and computational time. I also included Carl's comments and created a few vectors outside the loop. 

library(Matrix)

medianpricemat <-  as.vector(medianpricemat)
D1 <- Diagonal(x=pricedraws)

elasticitylist = foreach(i=1:nhousehold) %dopar% {
    pricedraws = out$betadraw[i,12,] 
    tmp = probarray[,,i] %*% D1
    elasticitydraws[,,i] = as.matrix(tmp) * medianpricemat
    MAPelasticity[i,] = rowMeans(elasticitydraws[,,i])
}

A less obvious hack is to avoid creation of the diagonal matrix and the matrix multiplication:

D2 <- rep(pricedraws, each=m)

elasticitylist = foreach(i=1:nhousehold) %dopar% {
    pricedraws = out$betadraw[i,12,] 
    tmp = probarray[,,i] * D2  # element wise multiplication 
    elasticitydraws[,,i] = as.matrix(tmp) * medianpricemat
    MAPelasticity[i,] = rowMeans(elasticitydraws[,,i]) 
}

Upvotes: 2

Related Questions