user3745220
user3745220

Reputation: 41

matrix multiplication in R (incredibly slow)

I have the following piece of code:

Y.hat.tr <- array(0,c(nXtr,2))  
  for (i in 1:nXtr){
    #print(i)
    Y.hat.tr[i,2] <- ktr[,i]%*%solve(K + a*In)%*%Ytr
    #Y.hat.tr[i,2] <- ktr[,i]%*%chol2inv(chol((K + a*In)))%*%Ytr
  }
  Y.hat.tr[,1] <- Ytr

My problem is that nXtr =300, and ktr is a 300x300 matrix. This routine takes approx 30 seconds to run in R version 3.0.1. I have tried various approaches to reduce the run time, but to no avail.

Any ideas would be gratefully received. If any other information is required please let me know

I have now taken the solve(K + a*In)%*%Ytr out of the loop, which has helped, but I was hoping to somehow vectorise this piece of code. Having thought about this for a while, and also after looking through various posts I cannot see how this can be done?

Upvotes: 0

Views: 677

Answers (1)

Greg Snow
Greg Snow

Reputation: 49640

Maybe I am missing something (and without sample or simulated data to test on it is harder to check), but isn't your loop equivalent to:

Y.hat.tr[,2] <- t(ktr) %*% solve(K + a*In) %*% Ytr

?

Removing the loop altogether and using internal vectorized code may speed things up.

Also, you are using solve with 1 argument, often you can speed things by using solve with 2 arguments (fewer internal calculations), something like:

t(ktr) %*% solve( K + a*In, Ytr )

Your loop is of the type called embarrassingly parallel, which means that if you want to keep the loop and are working on a computer with more than 1 core (or have easy access to a cluster) then you could use the parallel package (and maybe simplest to convert using the foreach package) to run the calculations in parallel which sometimes can greatly speed up the process.

Upvotes: 1

Related Questions