0到n之间的所有k个数的组合，其和等于n，速度优化

Question

我有这个 R 函数来生成一个矩阵，该矩阵由 0 到 n 之间的 k 个数字的所有组合组成，其总和等于 n。这是我的程序的瓶颈之一，因为即使数量很少，它也会变得非常慢（因为它正在计算幂集）

这里是代码

sum.comb <-
function(n,k) {

 ls1 <- list()                           # generate empty list
 for(i in 1:k) {                        # how could this be done with apply?
    ls1[[i]] <- 0:n                      # fill with 0:n
 }
 allc <- as.matrix(expand.grid(ls1))     # generate all combinations, already using the built in function
 colnames(allc) <- NULL
 index <- (rowSums(allc) == n)       # make index with only the ones that sum to n
 allc[index, ,drop=F]                   # matrix with only the ones that sum to n
 }

Answer 1

除非你回答我关于n 和k 的典型值的问题，否则很难判断它是否有用（请这样做。）这是一个使用递归的版本，它似乎比 josilber 使用他的更快基准测试：

sum.comb3 <- function(n, k) {

   stopifnot(k > 0L)

   REC <- function(n, k) {
      if (k == 1L) list(n) else
      unlist(lapply(0:n, function(i)Map(c, i, REC(n - i, k - 1L))),
             recursive = FALSE)
   }

   matrix(unlist(REC(n, k)), ncol = k, byrow = TRUE)
}

microbenchmark(sum.comb(3, 10), sum.comb2(3, 10), sum.comb3(3, 10))
# Unit: milliseconds
#              expr      min       lq   median       uq      max neval
#  sum.comb2(3, 10) 39.55612 40.60798 41.91954 44.26756 70.44944   100
#  sum.comb3(3, 10) 25.86008 27.74415 28.37080 29.65567 34.18620   100

Answer 2

这是一种不同的方法，它在每次迭代时将集合的大小从 1 逐步扩展为 k，修剪总和超过 n 的组合。如果 k 相对于 n 较大，这应该会导致加速，因为您不需要计算接近幂集大小的任何内容。

sum.comb2 <- function(n, k) {
  combos <- 0:n
  sums <- 0:n
  for (width in 2:k) {
    combos <- apply(expand.grid(combos, 0:n), 1, paste, collapse=" ")
    sums <- apply(expand.grid(sums, 0:n), 1, sum)
    if (width == k) {
      return(combos[sums == n])
    } else {
      combos <- combos[sums <= n]
      sums <- sums[sums <= n]
    }
  }
}

# Simple test
sum.comb2(3, 2)
# [1] "3 0" "2 1" "1 2" "0 3"

这是一个小 n 和大 k 的加速示例：

library(microbenchmark)
microbenchmark(sum.comb2(1, 100))
# Unit: milliseconds
#               expr      min      lq   median       uq      max neval
#  sum.comb2(1, 100) 149.0392 158.716 162.1919 174.0482 236.2095   100

这种方法在一秒钟内运行，而使用幂集的方法当然永远不会超过对 expand.grid 的调用，因为您最终会在结果矩阵中得到 2^100 行。

即使在不那么极端的情况下，在 n=3 和 k=10 的情况下，与原始帖子中的函数相比，我们也看到了 20 倍的加速：

microbenchmark(sum.comb(3, 10), sum.comb2(3, 10))
# Unit: milliseconds
#              expr       min        lq    median        uq       max neval
#   sum.comb(3, 10) 404.00895 439.94472 446.67452 461.24909 574.80426   100
#  sum.comb2(3, 10)  23.27445  24.53771  25.60409  26.97439  65.59576   100

Answer 3

请参阅partitions 包，尤其是compositions() 和blockparts()，它们作为整个矩阵生成器和迭代操作都会更快。如果这还不够快，请参阅关于组合和分区生成算法（无环、格雷码和并行）的各种出版物，例如 Daniel Page's research。

library(partitions)
library(microbenchmark)

# rcpp_comps is an Rcpp implementation of compositions using loop
# free grey code, just for illustrative purposes.

# Just get the matrix
microbenchmark( compositions(3,10), compositions(10,3),
                blockparts(rep(10,3),10), blockparts(rep(3,10),3),
                rcpp_comps(10), times=10)
## Unit: microseconds
##                        expr    min     lq median   uq    max neval
##         compositions(3, 10) 1967.4 2050.9 2097.1 2173 3189.6    10
##         compositions(10, 3)  618.2  638.5  654.6  688  700.7    10
##  blockparts(rep(10, 3), 10)  612.2  620.8  645.6  663  963.5    10
##   blockparts(rep(3, 10), 3) 2057.2 2089.2 2176.0 2242 3116.4    10
##              rcpp_comps(10)  359.9  360.7  367.6  378  404.2    10

Answer 4

使用 lapply 可以完成以下操作。

ls1 <- list()
for(i in 1:k) {
  ls1[[i]] <- 0:n
}

试着用这个代替，看看你是否能加快速度。

ls1 = lapply(1:k,function(x) 0:n)

我将 'ls' 更改为 'ls1' 因为 ls() 是一个 R 函数。

Answer 5

那么短的东西呢：

comb = function(n, k) {
    all = combn(0:n, k)
    sums = colSums(all)
    all[, sums == n]
}

然后是这样的：

comb(5, 3)

根据您的要求生成矩阵：

     [,1] [,2]
[1,]    0    0
[2,]    1    2
[3,]    4    3

---

感谢 @josilber 和原始发帖人指出 OP 需要所有重复排列而不是组合。排列的类似方法如下所示：

perm = function(n, k) {
    grid = matrix(rep(0:n, k), n + 1, k)
    all = expand.grid(data.frame(grid))
    sums = rowSums(all)
    all[sums == n,]
}

然后是这样的：

perm(5, 3)

根据您的要求生成一个矩阵：

    X1 X2 X3
6    5  0  0
11   4  1  0
16   3  2  0
21   2  3  0
26   1  4  0
31   0  5  0
...