R：查找所有组合而不替换稀疏矩阵答案

【问题标题】：R: Find all combinations without replacement of a sparse matrixR：查找所有组合而不替换稀疏矩阵
【发布时间】：2019-09-01 03:18:05
【问题描述】：

我想找到一个大稀疏矩阵的所有可能组合（无需替换）。每个组合最多只能从每一行和每一列中选择一次。我的目标是找到使所选条目的总和最大化的组合。

假设我有以下矩阵：

6 8 . .
. 5 7 .
. 6 . 9

有 4 种可能的组合（就 i 和 j 而言）：[(1,1),(2,2),(3,4)],[(1,1),(2,3), (3,2)],[(1,2),(2,3),(3,2)],[(1,2),(2,3),(3,4)]

我的结果应该是每个可能组合的条目总和，我的最终目标是找到使该结果最大化的组合 ([(1,2),(2,3),(3,4)] =在本例中为 8 + 7 + 9 = 24）。

编辑：这里是生成稀疏矩阵的完整代码，我想找到它的最佳组合：

library(data.table)
library(ggplot2)
library(haven)
library(Matrix)
library(evd) 

set.seed(12345)

N1 <- 100
M <- 100
I1 <- 10
I2 <- 2
I <- I1 * I2
N <- N1 * I2
J <- 5
p_c_A = 0.02
p_c_B = 0.01
p_0 = 0.05
p_1 = 0.2

dt_workers<- data.table(worker_id = 1:N, 
                           firm_id = sample.int(M, N, replace = TRUE),
                           worker_type = sample.int(I1, N, replace = TRUE)) 

dt_workers[, worker_ethnicity := 1 * (worker_id > N1)]

dt_firms <- data.table(firm_id = 1:M, 
                         firm_type = sample(J) )


sys_util <- matrix(NA, nrow=I1, ncol=J)
for(i in 1:dim(sys_util)[1]){
  for(j in 1:dim(sys_util)[2]){
    sys_util[i,j] <- i * j}
}


joint_surplus

con_A <- matrix(data = runif(N1 * M), nrow = N1, ncol = M)
con_B <- matrix(data = runif(N1 * M), nrow = N1, ncol = M)

con_A <- 1 * (con_A < p_c_A)
con_B <- 1 * (con_B < p_c_B)

p_meet_A <- con_A * p_1 + (1 - con_A) * p_0
p_meet_B <- con_B * p_1 + (1 - con_B) * p_0

meet_A <- matrix(data = runif(N1 * M), nrow = N1, ncol = M)
meet_B <- matrix(data = runif(N1 * M), nrow = N1, ncol = M)

meet_A <- 1* ( meet_A < p_meet_A )
meet_B <- 1* ( meet_B < p_meet_B )

meet <- rbind(meet_A,meet_B)

meet_sparse <- Matrix(meet, sparse = TRUE)
util <- which (meet_sparse>0, arr.ind=T)

n_draws <- dim(util)[1]


mu = 0
sigma = 10
idio = rgumbel(n=n_draws, loc=mu, scale=sigma)

util <- cbind(util,idio)
sys <- vector()
for(k in 1:dim(util)[1]){
  g <- util[k,1]
  f <- util[k,2]
  i <- dt_workers[g, worker_type]
  j <- dt_firms[f, firm_type]
  sys[k] = sys_util[i,j]
}
util <- cbind(util,sys)

total_util = util[,3] + util[,4]

M <- sparseMatrix(
  i = util[,1],
  j = util[,2],
  x = total_util
)
dat <- as.data.frame(summary(M))
dat <-dat[order(dat$i, dat$j),]
rownames(dat) <- NULL

【问题讨论】：

提供一些示例输入输出
您只想要返回最大值、采样的矩阵位置还是返回的采样值？
我想要最大化值总和的矩阵位置。
思考：我认为它应该被表述为一个线性规划问题
是的，你可能是对的。谢谢！

标签： r combinations sparse-matrix combinatorics

【解决方案1】：

library(Matrix)

M <- sparseMatrix(
  i = c(1, 1, 2, 2, 3, 3),
  j = c(1, 2, 2, 3, 2, 4),
  x = c(6, 8, 5, 7, 6, 9)
)

#> 3 x 4 sparse Matrix of class "dgCMatrix"
#>             
#> [1,] 6 8 . .
#> [2,] . 5 7 .
#> [3,] . 6 . 9

dat <- as.data.frame(summary(M))

#>   i j x
#> 1 1 1 6
#> 2 1 2 8
#> 3 2 2 5
#> 4 3 2 6
#> 5 2 3 7
#> 6 3 4 9

row_indices <- unique(dat$i)
col_indices <- split(dat$j, dat$i)

#> $`1`
#> [1] 1 2
#> 
#> $`2`
#> [1] 2 3
#> 
#> $`3`
#> [1] 2 4

all_combinations_with_atmost_one_per_row <- do.call(expand.grid, col_indices) 

#>   1 2 3
#> 1 1 2 2
#> 2 2 2 2
#> 3 1 3 2
#> 4 2 3 2
#> 5 1 2 4
#> 6 2 2 4
#> 7 1 3 4
#> 8 2 3 4

more_than_one_per_col <- apply(all_combinations_with_atmost_one_per_row, MARGIN = 1, anyDuplicated)

#> [1] 3 2 0 3 0 2 0 0

combinations <- all_combinations_with_atmost_one_per_row[!more_than_one_per_col, , drop = FALSE]

#>   1 2 3
#> 3 1 3 2
#> 5 1 2 4
#> 7 1 3 4
#> 8 2 3 4

lapply(
  split(combinations, 1:nrow(combinations)),
  function(cols) {
    elements <- data.frame(i = row_indices, j = unlist(cols))
    elements$value <- M[as.matrix(elements)]
    list(elements = elements, sum = sum(elements$value))
  }
)
#> $`1`
#> $`1`$elements
#>   i j value
#> 1 1 1     6
#> 2 2 3     7
#> 3 3 2     6
#> 
#> $`1`$sum
#> [1] 19
#> 
#> 
#> $`2`
#> $`2`$elements
#>   i j value
#> 1 1 1     6
#> 2 2 2     5
#> 3 3 4     9
#> 
#> $`2`$sum
#> [1] 20
#> 
#> 
#> $`3`
#> $`3`$elements
#>   i j value
#> 1 1 1     6
#> 2 2 3     7
#> 3 3 4     9
#> 
#> $`3`$sum
#> [1] 22
#> 
#> 
#> $`4`
#> $`4`$elements
#>   i j value
#> 1 1 2     8
#> 2 2 3     7
#> 3 3 4     9
#> 
#> $`4`$sum
#> [1] 24

^{由reprex package (v0.2.1) 于 2019-04-10 创建}

最好的组合是res[[which.max(sapply(res, `[[`, "sum"))]]

$elements
  i j value
1 1 2     8
2 2 3     7
3 3 4     9

$sum
[1] 24

【讨论】：

非常感谢您的详细解答！出于某种原因，您的代码在我的模拟中不起作用。问题从命令 col_indices
我想我在代码中发现了问题。它应该是： col_indices
它适用于非常小的数据集。当数据有点大时，没有足够的内存来执行该步骤：all_combinations_with_atmost_one_per_row
我的回答不管用。请参阅M <- sparseMatrix( i = c(1, 2, 2), j = c(1, 1, 2), x = c(1, 3, 1) )，其中最佳解决方案是[(2, 1)]，总和为3，而我发现[(1, 1), (2, 2)] 的总和为2。
是的，您需要添加值为 0 的不匹配选项（在您的示例中为 [1, 0]）（婚姻问题术语中的“单身”）。我想将选项 i=0 和 j=0 添加到基本矩阵很容易，但没有限制该原始/列之和应小于 1。

【解决方案2】：

我在 Aurèle 的帮助下找到了使用线性规划的解决方案：

f.con <- matrix(,nrow = dim(dat)[1],ncol=0)
for(k in 1: N){   
    vec <- 1 * (dat[,1] == k)
    f.con <- cbind(f.con , vec )
}

for(k in 1: M){   
  vec <- 1 * (dat[,2] == k)
  f.con <- cbind(f.con , vec )
}
f.con <- t(f.con)
f.obj <- dat[,3]
f.dir <- rep ("<=", dim(f.con)[1])
f.rhs <- rep (1, dim(f.con)[1])
res = lp (direction = "max", f.obj, f.con, f.dir, f.rhs ,  all.int=TRUE)$solution

【讨论】：