有效地选择整数的组合

Question

假设我们有一个用 0 填充的 5x5 矩阵。

myMatrix <- matrix(rep(0, 25), ncol = 5)

现在，让我们选择 1 到 5 之间的三元组整数。

triplet <- c(1,2,3)

对于这个三元组的所有组合，我们现在在矩阵中加 1，使用此功能：

addCombinationsToMatrix <- function(.matrix, .triplet){
    indexesToChange <- as.matrix(expand.grid(.triplet, .triplet))
    .matrix[indexesToChange] <- .matrix[indexesToChange] + 1
    .matrix
}

使用函数，我们从

myMatrix

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

至

myMatrix <- addCombinationsToMatrix(myMatrix, triplet)
myMatrix

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

如果我们选择另一个三胞胎，我们将继续

nextTriplet <- 2:4
myMatrix <- addCombinationsToMatrix(myMatrix, nextTriplet)
myMatrix

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

因此，行列组合表示两个整数出现的频率三联体：3 和 4 一起出现过一次，2 和 3 出现过一次一起出现过两次。

问题: 如何挑选三胞胎，使得每个组合（1-2、1-3、1-4 ...）至少被选中一次并且三胞胎的数量被最小化。

我正在寻找一种算法来选择下一个三元组。

理想情况下可以扩展到

任意大矩阵（10x10、100x100 ...）
任意大向量（四元组、五元组、n-元组）
必须至少选择组合的任意次数

示例：

myMatrix
myMatrix <- addCombinationsToMatrix(myMatrix, 1:3)
myMatrix
myMatrix <- addCombinationsToMatrix(myMatrix, 3:5)
myMatrix
myMatrix <- addCombinationsToMatrix(myMatrix, c(1,4,5))
myMatrix
myMatrix <- addCombinationsToMatrix(myMatrix, c(2,4,5))
myMatrix

编辑：可以肯定的是：答案不一定是 R 代码。它也可以是其他语言，甚至是伪代码。

编辑 2：我现在想到，有不同的方法来衡量效率。我的意思是，该算法应该尽可能少地进行迭代。算法快也很酷，但不是这里的主要目标。

Answer 1

这是一个使用 data.table 跟踪矩阵计数并使用 RcppAlgos 生成组合的选项：

library(RcppAlgos)
library(data.table)

M <- 100 #5 #10 #100
sz <- 5 #3 #4 5 
minpick <- 3 #1 #2
d <- integer(M)

system.time({
    universe <- as.data.table(comboGeneral(M, 2L, nThreads=4L))[, count := 0L]
    ntuples <- 0
    while (universe[, any(count < minpick)]) {
        v <- universe[order(count), head(unique(c(V1[1L:2L], V2[1L:2L])), sz)]
        universe[as.data.table(comboGeneral(v, 2L, nThreads=4L)), on=.NATURAL, count := count + 1L]
        ntuples = ntuples + 1L
    }
    ntuples
})
#   user  system elapsed 
#  26.82    9.81   28.75 

m <- matrix(0L, nrow=M, ncol=M)
m[as.matrix(universe[, V1:V2])] <- universe$count
m + t(m) + diag(d)

这是一个贪心算法，因此，我不确定这是否会产生最少数量的元组。

Answer 2

好问题！这出现在调查设计中，您需要几个不同版本的调查，每个版本只包含问题的一个子集，但您希望每对（或 t 元组）问题至少被问过一次。

这称为覆盖设计，是经典集合覆盖问题的变体。正如您在优秀的 Mathematics Stack Exchange post on the topic, folks use notation C(v, k, t) indicating the minimum number of k-element subsets you need to draw (k=3 in your case) from a v-element set (v=5 in your case) such that every t-element subset in the entire set (t=2 in your case) is contained within one of your selected subsets. Folks have evaluated this function for many different (v, k, t) tuples; see, for instance, https://ljcr.dmgordon.org/cover/table.html 中所读到的那样。我们可以从 table 中读出 C(5, 3, 2) = 4，以下是一种可能的设计：

首先，这个问题是 NP-hard 问题，因此所有已知的精确算法都将在输入 v、k 和 t 中呈指数缩放。因此，虽然您可以通过枚举或一些更聪明的精确方法（例如整数规划）精确地解决小实例，但随着问题规模变得非常大，您可能需要求助于启发式方法。

这个方向的一种可能性是字典覆盖，如 https://arxiv.org/pdf/math/9502238.pdf 中所提议的（您会注意到上面链接的站点上的许多解决方案列出了 "lex covering" 作为构造方法）。基本上，您按字典顺序列出所有可能的 k 元组：

然后你贪婪地添加覆盖之前发现最多的 t 元组的 k 元组，使用字典顺序打破平局。

算法在我们的例子中是这样工作的：

开始时每个三元组包含三个不同的二元组，所以我们添加 123 因为它在字典上最早。
这样就覆盖了12、13、23的二元组，剩下的二元组都没有覆盖.许多 3 元组覆盖了 3 个以上的 2 元组，例如145 和 245。我们选择 145，因为它在字典序上排在第一位，涵盖 14、45 和 15。
现在我们还有 4 个未覆盖的二元组 -- 24、25、34 和 35。没有 3 元组覆盖其中的 3 个，但有几个覆盖 2 个，例如234 和 345。我们 select 234 作为字典编年史最早的
我们还有两个未覆盖的二元组 -- 25 和 35。我们 select 235 作为唯一涵盖两者的三元组。

我们最终得到了如上所示的精确解。重要的是，这只是一种启发式方法——它不能保证 4 是覆盖具有 5 个元素的集合中所有对所需的最小三元组数。在这种情况下，Schönheim 的下界（在上面的链接文章中提供了参考）使我们相信，实际上 C(5, 3, 2) 不能小于 4。我们得出结论，字典序覆盖的解是实际上是最优的。

您需要进行调整以将每个 t 元组覆盖一定次数 r。一个明显的方法就是重复每个要覆盖的元组 "r" 次，然后像往常一样运行 lex 覆盖（例如，在上面的第一步中，每个 3 元组将覆盖 9 个 2 元组r=3).当然，由于使用了 lex 覆盖，这仍然是您整体问题的启发式方法。

Answer 3

由于这个问题要求的是覆盖设计的算法方法，我将提供一种使用 R 中的整数规划给出准确答案（也称为最佳设计）的方法。对于您正在考虑的每个 k 元组（k =3 对你来说，因为你选择的是三元组），定义一个决策变量，如果你将它包含在你的设计中则取值为 1，否则取值为 0。因此，在您的情况下，您将定义 x_123 以指示是否选择了元组 (1,2,3)，x_345 表示 (3,4,5)，依此类推。

优化模型的 objective 是最小化所选元组的数量，也就是所有决策变量的总和。但是，对于每个 t 元组（在您的情况下为 t=2），您需要包含一个包含该 t 元组的决策变量。这会为每个 t 元组产生一个约束。例如，我们将 x_123+x_124+x_125 >= 1 作为约束，要求 12 对位于某个选定的元组中。

这会产生以下优化模型：

min  x_123+x_124+...+x_345
s.t. x_123+x_124+x_125 >= 1  # constraint for 12
     x_123+x_134+x_135 >= 1  # constraint for 13
     ...
     x_145+x_245+x_345 >= 1  # constraint for 45
     x_ijk binary for all i, j, k

您可以将其扩展为要求每个 t 元组重复 r 次，方法是将每个不等式的右侧更改为 "r" 并要求所有变量都是整数而不是二进制。

这很容易用 R 中的 lpSolve 这样的包解决：

library(lpSolve)
C <- function(v, k, tt, r) {
  k.tuples <- combn(v, k)
  t.tuples <- combn(v, tt)
  mod <- lp(direction="min",
            objective.in=rep(1, ncol(k.tuples)),
            const.mat=t(apply(t.tuples, 2, function(x) {
              apply(k.tuples, 2, function(y) as.numeric(sum(x %in% y) == tt))
            })),
            const.dir=rep(">=", ncol(t.tuples)),
            const.rhs=rep(r, ncol(t.tuples)),
            all.int=TRUE)
  k.tuples[,rep(seq_len(ncol(k.tuples)), round(mod$solution))]
}
C(5, 3, 2, 1)
#      [,1] [,2] [,3] [,4]
# [1,]    1    1    1    3
# [2,]    2    2    2    4
# [3,]    3    4    5    5
C(5, 3, 2, 3)
#      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,]    1    1    1    1    1    1    2    2    2     3
# [2,]    2    2    2    3    3    4    3    3    4     4
# [3,]    3    4    5    4    5    5    4    5    5     5

虽然这可以准确解决您的问题，但它不能很好地扩展到大问题规模。这是因为问题是 NP 难的——没有已知的精确算法可以很好地扩展。如果您需要解决大型问题实例，那么此处其他答案中推荐的启发式方法是您的最佳选择。或者你可以用整数规划来解决（就像我们在这里做的那样）并设置超时；那么您将使用超时找到的最佳解决方案，这是对整个问题的启发式解决方案。

有效地选择整数的组合

Efficiently picking combinations of Integers

algorithm

r

mathematical-optimization