【问题标题】：Multiplication of many matrices in RR中许多矩阵的乘法
【发布时间】：2019-04-16 13:50:51
【问题描述】：

我想将几个相同大小的矩阵与一个初始向量相乘。在下面的示例中，p.state 是 m 元素的向量，tran.mat 是列表，其中每个成员都是 m x m 矩阵。

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

上面的代码给出了正确的答案，但是当length(tran.mat) 很大时可能会很慢。我想知道是否有更有效的方法来做到这一点？

下面是一个带有m=3 和length(mat)=10 的示例，可以生成这个：

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)

注意：tran.mat 不必是一个列表，它只是当前写成一个。

在几个 cmets 后编辑：

Reduce 在m 很小的时候很有用。但是，当m=6 循环执行上述两种解决方案时。库（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)
  }

)

这会导致

        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

【问题讨论】：

标签： r loops matrix multiplication markov-chains

【解决方案1】：

更快的方法是使用Reduce() 对矩阵列表进行顺序矩阵乘法。

通过这种方式，您可以获得大约 4 倍的加速。下面是一个测试代码示例，列表中有 1000 个元素而不是 10 个元素，以便更轻松地查看性能改进。

代码

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)
  }
)

输出

    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

你可以看到 reduce 方法快了大约 3.8 倍。

【讨论】：

PSA：不要用字符串引用 R 名称，用反引号转义它们（即不要写 '%*%'，写 `%*%`）。 R 甚至允许使用字符串引用的函数名称这一事实是一种颠覆性的类型系统。除此之外，为什么不使用Reduce(`%*%`, tran.mat, p.state)？
当m 很大时，上述方法不起作用。查看原始问题中的编辑
我怀疑这是因为乘法的顺序。即在循环中，我们将1 x m 矩阵乘以m x m 矩阵length(tran.mat) 次，而使用Reduce 我们将m x m 矩阵乘以length(tran.mat) 次。前者的操作要少得多。

【解决方案2】：

我不确定这会更快，但它会更短：

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

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

注意

我们使用了这个测试输入：

m <- matrix(1:9, 3)
L <- list(m^0, m, m^2, m^3)

【讨论】：

你也是，在字符串中引用函数名？ :(
我强烈反对。它颠覆了类型系统并产生了语义后果。句法差异（或者，更确切地说，相似性）纯粹是偶然的，事实上，相当具有误导性。您不会“调用字符串”，单纯的概念是无稽之谈。而 R 在幕后发挥作用（通过match.fun 和语法松散）这一事实对初学者和中级 R 用户造成了极大的伤害，因为他们误解了正在发生的事情。根据您的论点，写"a" = 1 是“好的”。你肯定同意那个代码是荒谬的。
这是一种逃避。我的“意见”是有道理的，基于事实，而不是一时兴起。这就是为什么我很惊讶看到像你这样经验丰富的 R 用户不同意的原因。将 all 函数名作为参数传递时是否引用它们？还是只有运营商？

【解决方案3】：

我将使用 package Rfast 中的一个函数来减少乘法的执行时间。不幸的是，for循环的时间不能减少。

名为Rfast::eachcol.apply 的函数是一个很好的解决方案。您的乘法也是函数crossprod，但对于我们的目的来说它很慢。

这里有一些辅助函数：

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
}

这是一个例子：

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

oper 参数用于对每个元素的操作，apply 用于对每列的操作。在大型矩阵中应该很快。我无法在我的笔记本电脑上测试它以获得更大的矩阵。

【讨论】：