寻找更快的 R 计算以进行迭代矩阵计算

Question

我有以下代码位于优化例程中。因此，虽然速度相当快，但分析显示产生结果的行被称为res我的代码中最大的瓶颈。

我尝试了很多方法来改进这一点，并最终得到了最后一行：

res <- exp(X %*% log(pr.t) + mX %*% log(1 - pr.t)) %*% wts

在我的问题中，矩阵的元素X是固定的，并且不会随着迭代而改变。因此，我还可以计算、存储和回收X和mX。每次迭代发生的变化是我在对象中计算的一些概率pr.t。

我尝试过 Rcpp，但 Rcpp 与我工作中的 R 代码一样快。

我现在向这个小组发出呼吁，看看是否有人能找到一种绝妙的方法来加快最终产品的生产线速度res。下面是设置问题的示例代码，给出了实际问题的可重现示例。

X <- matrix(sample(c(0,1), 5000, replace = TRUE), 1000, 5)
mX <- 1 - X
pr.t <- matrix(runif(75), 5, 15)
wts <- runif(15)
res <- exp(X %*% log(pr.t) + mX %*% log(1 - pr.t)) %*% wts

Answer 1

首先，如果nrow(X)远大于nrow(unique(X))，如您的示例所示，您只需计算nrow(unique(X))的可能值res，然后计算索引。这将减少乘法的大小。其次，exp、log和前两个矩阵乘法可以替换为直接执行乘法的函数：

n <- ncol(X)
X0 <- unique(X)
mX0 <- 1 - X0 # only needed for `res2`
i <- match(X %*% 2^((n - 1):0), X0 %*% 2^((n - 1):0))
X02 <- collapse::setop(X0*n, "+", 1:n, rowwise = TRUE)
mode(X02) <- "integer"

f <- function(pr.t, wts) {
  pr.t2 <- rbind(1 - pr.t, pr.t)
  res <- pr.t2[X02[,1],]
  if (n != 1L) for (j in 2:n) res <- res*pr.t2[X02[,j],]
  (res %*% wts)[i,,drop = FALSE]
}

microbenchmark::microbenchmark(
  res = exp(X %*% log(pr.t) + mX %*% log(1 - pr.t)) %*% wts,
  res2 = (exp(X0 %*% log(pr.t) + mX0 %*% log1p(-pr.t)) %*% wts)[i,,drop = FALSE],
  res3 = f(pr.t, wts),
  check = "equal",
  unit = "relative"
)
#> Unit: relative
#>  expr       min        lq      mean    median        uq       max neval
#>   res 47.141593 44.735632 27.885564 40.779310 33.565934 7.1477733   100
#>  res2  2.088496  1.900383  1.317103  1.751724  1.483516 0.8380567   100
#>  res3  1.000000  1.000000  1.000000  1.000000  1.000000 1.0000000   100