R组合: 探寻比基本R更快更高效的方法（包括软件包、代码和并行CPU）。

Question

我正在使用基本的 R 来进行组合。

例如，假设我有一个 2 行 5 列的矩阵：

 z<-matrix(c(1, 2, 1, 3, 2, 2, 1, 3, 2, 1),nrow=2,ncol=5,byrow = TRUE)

[,1] [,2] [,3] [,4] [,5]

[1,]    1    2    1    3    2

[2,]    2    1    3    2    1

我正在使用下面的代码来生成从5列中选3列的组合：

l<- apply(X = combn(seq_len(ncol(z)), 3),MAR = 2,FUN = function(jj) {apply(z[, jj], 1, paste, collapse="") })

这将导出我所需的内容：

[,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]  [,10]

[1,] "121" "123" "122" "113" "112" "132" "213" "212" "232" "132"

[2,] "213" "212" "211" "232" "231" "221" "132" "131" "121" "321"

问题出现在我使用矩阵中的大数据时，例如当我有一个包含15000行和17列的矩阵，并且我需要从这17列中获取10个集合的组合时。在这个例子中，导出需要很长时间。

对于这个组合的例子，是否有比基本的R更快更有效的方法（也许是一些软件包或代码，或者使用并行CPU）？

我正在使用Windows 7 64位，FX 8320，16GB RAM。

- pshls

---

不知道你能节省多少时间，但你可以简化一下你的代码：apply(z,1,function(x) combn(x,3,FUN=paste,collapse="")) 将会产生 t(l)。 - nicola

我有一种感觉，我们会看到Dirk过来推荐Rcpp。 :) 如果你还没有尝试过，那可能是一个不错的选择去探索。 - Alex A.

对于您的示例，您希望生成292百万个组合（17选10乘以15000），所以花费一些时间并不令人意外... - josliber

@nicola 我建议您将此作为答案添加，并附上一些基准测试数据 - 我发现在选择n列的100 x 17矩阵时，您的代码运行时间为0.3秒，而OP的代码则需要16秒。 - josliber

当我处理生成组合时，我注意到paste函数会减慢代码的速度。将数据保留在矩阵形式中可以使代码运行更加高效。 - inscaven

Answer 1

正如@inscaven所指出的，真正的时间压力来自于paste。如果我们只需要生成所有17个选10个组合15000次，那么使用几个高度优化的包，如R、arrangements和RcppAlgos(我是作者)，这并不需要太长时间。

set.seed(101)
testMat <- matrix(sample(1000, 15000 * 17, TRUE), nrow = 15000)

library(arrangements)
system.time(lapply(1:15000, function(x) {
    temp <- combinations(x = testMat[x, ], k = 10)
    x
}))
  user  system elapsed 
 6.879   2.133   9.014

library(RcppAlgos)
system.time(lapply(1:15000, function(x) {
    temp <- comboGeneral(testMat[x, ], 10)
    x
}))
  user  system elapsed 
 5.770   2.178   7.953

与在基础R中加载的combn相比：

system.time(lapply(1:15000, function(x) {
    temp <- combn(testMat[x, ], 10)
    x
}))
    user  system elapsed 
 261.163   1.093 262.608 

如果我们必须将结果合并成字符矩阵，那么在基本的R语言中我们几乎无法做更多的工作。即使使用上述任何一种优化库，我们仍然需要循环遍历所有行并将结果粘贴在一起，这会导致速度变慢。

system.time(t1 <- lapply(1:50, function(x) {
    combn(testMat[x, ], 10, paste0, collapse = "")
}))
  user  system elapsed 
 6.847   0.070   6.933

## from package arrangements
system.time(t2 <- lapply(1:50, function(x) {
    apply(combinations(x = testMat[x, ], k = 10), 1, paste0, collapse = "")
}))
  user  system elapsed 
 6.318   0.032   6.353

这并不是真正的胜利。我们需要一种新的方法。 引入Rcpp

//[[Rcpp::export]]
CharacterVector pasteCombos(int n, int r, CharacterVector v, int numRows) {

    int r1 = r - 1, r2 = r - 2;
    int numIter, count = 0;
    CharacterVector comboVec = Rcpp::no_init_vector(numRows);

    std::vector<int> z(r);
    std::iota(z.begin(), z.end(), 0);

    while (count < numRows) {
        numIter = n - z[r1];
        if ((numIter + count) > numRows)
            numIter = numRows - count;

        for (int i = 0; i < numIter; ++i, ++count, ++z[r1])
            for (int k = 0; k < r; ++k)
                comboVec[count] += v[z[k]];

        for (int i = r2; i >= 0; i--) {
            if (z[i] != (n - r + i)) {
                ++z[i];
                for (int k = (i + 1); k < r; ++k) 
                    z[k] = z[k - 1] + 1;

                break;
            }
        }
    }

    return comboVec;
}

该函数简单地生成所有的v个选择r个的组合，并通过+=实时粘贴结果。这样就可以生成一个无需处理矩阵行的向量。让我们看看是否有任何改进。

numCombs <- choose(17, 10)
charMat <- matrix(as.character(testMat), nrow = 15000)

funOP <- function(z, r) {
    apply(X = combn(seq_len(ncol(z)), r), MAR = 2,FUN = function(jj) {apply(z[, jj], 1, paste, collapse="") })
}

system.time(t1 <- funOP(testMat[1:100, ], 10))
   user  system elapsed 
 22.221   0.110  22.330 

system.time(t2 <- lapply(1:100, function(x) {
     pasteCombos(17, 10, charMat[x,], numCombs)
}))
  user  system elapsed 
 7.890   0.085   7.975

近乎3倍的速度提升...不错，但我们可以做得更好。

进入parallel

library(parallel)
system.time(t3 <- mclapply(1:100, function(x) {
    pasteCombos(17, 10, charMat[x,], numCombs)
}, mc.cores = 8)) ## you will have to adjust this on your computer.. I'm running MacOS with 8 cores
  user  system elapsed 
 1.430   0.454   1.912

现在我们才开始！快了近12倍！！

这里进行一次健康检查：

all.equal(t1, do.call(rbind, t2))
# [1] TRUE
all.equal(t1, do.call(rbind, t3))
# [1] TRUE

假设我们可以在2秒内完成100行，那么总共完成任务需要 2 * 150 = 300 秒 = 5 分钟。

- Joseph Wood