R: Найти все комбинации без замены разреженной матрицы

Question

Я хочу найти все возможные комбинации (без замены) большой разреженной матрицы. Каждую комбинацию можно выбрать не более одного раза из каждой строки и столбца. Моя цель - найти комбинацию, которая максимизирует сумму выбранных записей.

Скажите, у меня есть следующая матрица:

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

Answer 1

Я нашел решение с помощью линейного программирования с помощью 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

Answer 2

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

^{Создано в 2019-04-10 пакетом Представления (v0.2.1)}

И лучшая комбинация найдена с 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