4

I want to multiply several matrices of the same size with an inital vector. In the example below p.state is vector of m elements and tran.mat is list where each member is an m x m matrix. 

for (i in 1:length(tran.mat)){
  p.state <- p.state %*% tran.mat[[i]]
}

The code above gives the correct answer but can be slow when length(tran.mat) is large. I was wondering if there was a more efficient way of doing this?

Below is an example with a m=3 and length(mat)=10 that can generate this:

p.state <- c(1,0,0)
tran.mat<-lapply(1:10,function(y){apply(matrix(runif(9),3,3),1,function(x){x/sum(x)})})

for (i in 1:length(tran.mat)){
  p.state <- p.state %*% tran.mat[[i]]
}

print(p.state)

NB: tran.mat does not have to be a list it is just currently written as one.

Edit after a few comments:

Reduce is useful when m is small. However when m=6 the loop out performed both the above solutions. library(rbenchmark)

p.state1 <- p.state <- c(1,0,0,0,0,0)
tran.mat<-lapply(1:10000,function(y){t(apply(matrix(runif(36),6,6),1,function(x){x/sum(x)}))})

tst<-do.call(c, list(list(p.state), tran.mat))

benchmark(
  'loop' = {
    for (i in 1:length(tran.mat)){
      p.state <- p.state %*% tran.mat[[i]]
    }
  },
  'reduce' = {
    p.state1 %*% Reduce('%*%', tran.mat)
   },
  'reorder' = {
    Reduce(`%*%`,tran.mat,p.state1)
  }

)

This results in

        test replications elapsed relative user.self sys.self user.child     sys.child
   1    loop          100    0.87    1.000      0.87        0         NA        NA
   2  reduce          100    1.41    1.621      1.39        0         NA        NA
   3 reorder          100    1.00    1.149      1.00        0         NA        NA
wearetheboro
  • 41
  • 2

3 Answers3

2

A faster way is to use Reduce() to do sequential matrix multiplication on the list of matrices.

You can get approximately a 4x speedup that way. Below is an example of your code tested, with 1000 elements in the list instead of 10 to see the performance improvement more easily.

Code

library(rbenchmark)

p.state <- c(1,0,0)
tran.mat<-lapply(1:1000,function(y){apply(matrix(runif(9),3,3),1,function(x){x/sum(x)})})

benchmark(
  'loop' = {
    for (i in 1:length(tran.mat)){
      p.state <- p.state %*% tran.mat[[i]]
    }
  },
  'reduce' = {
    p.state %*% Reduce('%*%', tran.mat)
  }
)

Output

    test replications elapsed relative user.self sys.self user.child sys.child
1   loop          100    0.23    3.833      0.23        0         NA        NA
2 reduce          100    0.06    1.000      0.07        0         NA        NA

You can see the reduce method is about 3.8 times faster.

qdread
  • 3,389
  • 19
  • 36
  • 1
    PSA: don’t string-quote R names, backtick-escape them (i.e. don’t write `'%*%'`, write `\`%*%\``). The fact that R even allows string-quoted function names is a type system subverting abomination. Apart from that, why not use `Reduce(\`%*%\`, tran.mat, p.state)`? – Konrad Rudolph Nov 13 '18 at 17:10
  • The above does not work when `m` is large. See the edit in the original question – wearetheboro Nov 14 '18 at 10:25
  • I suspect that this is because of the order of the multiplication. i.e in the loop we're multiplying an `1 x m` matrix by an `m x m` matrix `length(tran.mat)` times whereas using `Reduce` we are multiplying `m x m` matricies `length(tran.mat)` times. The former has a lot less operations. – wearetheboro Nov 14 '18 at 10:30
1

I am not sure that this will be any faster but it is shorter:

prod <- Reduce("%*%", L)

all.equal(prod, L[[1]] %*% L[[2]] %*% L[[3]] %*% L[[4]])
## [1] TRUE

Note

We used this test input:

m <- matrix(1:9, 3)
L <- list(m^0, m, m^2, m^3)
G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341
  • You too, quoting function names in strings? :( – Konrad Rudolph Nov 13 '18 at 17:12
  • I strongly disagree. It subverts the type system and has semantic consequences. The syntactic difference (or, rather, similarity) is purely incidental and, in fact, quite misleading. You don’t “call a string”, the mere concept is nonsense. And the fact that R does magic (by way of `match.fun` and syntactic laxity) behind the scenes does beginners and intermediate R users a huge disservice because they misunderstand what’s going on. By your argument it’s “fine” to write `"a" = 1`. Surely you agree that that code is absurd. – Konrad Rudolph Nov 13 '18 at 17:26
  • That’s a cop-out. My “opinion” is reasoned and based on facts, not on a whim. Which is why I’m so surprised to see a seasoned R user such as yourself disagreeing. Do you quote *all* function names when passing them as arguments? Or only operators? – Konrad Rudolph Nov 13 '18 at 17:33
1

I am going to use a function from package Rfast to reduce the execution time of multiplication. Unfortunately, for loop's time can not be reduced.

The function called Rfast::eachcol.apply is a great solution for your purpose. Your multiplication is also the function crossprod but it is slow for our purpose.

Here are some helper functions:

mult.list<-function(x,y){
    for (xm in x){
        y <- y %*% xm
    }
    y
}

mult.list2<-function(x,y){
    for (xm in x){
        y <- Rfast::eachcol.apply(xm,y,oper="*",apply="sum")
    }
    y
}

Here is an example:

x<-list()
y<-rnomr(1000)
for(i in 1:100){
    x[[i]]<-Rfast::matrnorm(1000,1000)
}


microbenchmark::microbenchmark(R=a<-mult.list(x,y),Rfast=b<-mult.list2(x,y),times = 10)
 Unit: milliseconds
     expr        min         lq        mean     median         uq        max neval
        R 410.067525 532.176979 633.3700627 649.155826 699.721086 916.542414    10
    Rfast 239.987159 251.266488 352.1951486 276.382339 458.089342 741.340268    10

all.equal(as.numeric(a),as.numeric(b))
[1] TRUE

The argument oper is for the operation on each element and the apply for the operation on each column. In large matrices should be fast. I couldn't test it in my laptop for bigger matrices.

Manos Papadakis
  • 564
  • 5
  • 17