for循环比你想象的更高效
您的n (A,B) 对相乘问题并不等同于通常意义上的张量乘法，尽管 whuber 提出了一种非常巧妙的方法，通过将 B 作为块堆叠在稀疏矩阵。
您曾说过要避免使用 for 循环，但是当高效编程时，for-loop 方法实际上非常具有竞争力，我建议您重新考虑一下。
我将使用与 whuber 相同的符号，A 的维度为 n x k，B 的维度为 k x k x n，例如：
n <- 1e4
k <- 3
A <- array(rnorm(k*n),c(n,k))
B <- array(rnorm(k*k*n),c(k,k,n))
一个简单有效的for循环解决方案应该是这样的
justAForLoop <- function(A,B) {
n <- nrow(A)
for (i in 1:n) A[i,] <- A[i,] %*% B[,,i]
A
}
产生一个 n x k 矩阵的结果。
我修改了whuber的f3函数加载Matrix包，否则sparseMatrix函数不可用。我的f3 版本比原来的版本快得多，因为我在返回结果之前已经消除了最后一个矩阵转置。
通过此修改，它会向justAForLoop 返回相同的数值结果。
f3 <- function(a, b) {
require(Matrix)
n <- dim(b)[3]
k <- dim(b)[2]
i0 <- (1:n-1)*k+1
i <- rep(i0, each=k)
j <- 1:(k*n)
aa <- sparseMatrix(i, j, x=c(t(a)), dims=c(n*k, n*k))
bb <- matrix(aperm(b, c(1,3,2)), nrow=n*k)
(aa %*% bb)[i0, ]
}
现在我在新的 R 会话中重新运行 whuber 的模拟：
> k <- 3
> n <- 1e6
> a <- matrix(runif(k*n), ncol=k)
> b <- array(runif(k^2*n), dim=c(k,k,n))
>
> system.time(c1 <- f1(a,b))
user system elapsed
3.40 0.09 3.50
> system.time(c3 <- f3(a,b))
Loading required package: Matrix
user system elapsed
1.06 0.24 1.30
> system.time(c4 <- justAForLoop(a,b))
user system elapsed
1.27 0.00 1.26
for-loop 方法实际上是最快的，差一点。它比依赖sapply 的f1 快得多。 （我的机器是运行 R 3.6.0 的具有 32Gb RAM 的 Windows 10 PC）。
如果我再次运行所有三个方法，那么f3 将成为最快的，因为这一次 Matrix 包已经在搜索路径中并且不必重新加载：
> system.time(c1 <- f1(a,b))
user system elapsed
3.23 0.04 3.26
> system.time(c3 <- f3(a,b))
user system elapsed
0.33 0.20 0.53
> system.time(c4 <- justAForLoop(a,b))
user system elapsed
1.28 0.01 1.30
但是f3 使用的 RAM 比 for 循环多。在我的电脑上，我可以使用n=1e8 成功运行justAForLoop，而f1 和f3 都内存不足而失败。
总结
直接的for循环方法比sapply高效得多。
对于 n=10,000 矩阵乘法的问题，运行 for 循环既简单又高效，耗时
对于介于 1-1000 万之间的 n，whuber 的稀疏矩阵解决方案开始表现出色，尤其是在 Matrix 包已经加载的情况下。
for 循环使用三种方法中最少的 RAM。对于n 在我的具有 32Gb RAM 的 PC 上的 1 亿，只有 for 循环方法有效。