In one of my applications, I need to perform following task in R: svector <- array(0, dim=c(5, 100000, 360), dimnames=list(c('s1','s2','s3','s4','s5'), NULL, NULL))) tmatrix <- array(0, dim=c(100000, 360, 5, 5), dimnames=list(NULL, NULL, c('s1','s2','s3','s4','s5'), c('s1','s2','s3','s4','s5')))) #My algorithms compute all the elements in tmatrix (transition matrix among states from time t to t+1, for all entities represented by l index) #do transition for all l for t=1 to 359 for (l in 1:100000) { for(t in 1:359) { svector [, l,t+1] <- tmatrix[l,t,,] %*% svector [,l,t] } } The double loops make computation slow. I have been trying to see I can treat the svector and tmatrix as tensors, and use mul.tensor in tensorR to vectorize the computation, but so far I ways get one message or another indicating my incorrect usage. Can tensorR or any other package be used here to simply the calculation? If it can, would you kindly give some sample code Thank you in advance!