Максимизируйте диагональ матрицы, переставляя столбцы в R

Question

Я работаю с квадратной матрицей в R, мы можем назвать ее mat и хотел бы переставить столбцы (т.е. изменить их порядок), чтобы максимизировать сумму диагональных элементов. Я хочу сделать это с помощью методов линейного программирования, то есть опираясь на пакет оптимизации lpSolve. Решения для кода, конечно, приветствуются, но если это не удастся, любая помощь, формулирующая ее как задачу линейного программирования, будет оценена.

Мой вопрос похож на этот: Перестановка столбцов квадратной двусторонней таблицы сопряженности (матрица), чтобы максимизировать его диагональ . Тем не менее, в этом и других вопросах, которые я нашел в SO, считается достаточным go построчное максимизация диагонального элемента в этой строке. Проблема в том, что что-то вроде

mat2 <- mat[,max.col(mat, 'first')]

не будет работать для меня: у вас могут быть ситуации, когда строка имеет несколько равных максимумов или где (скажем) в строке X вы выбираете 11 по диагонали, а не 10, но, следовательно, в строке X + 1 вы должны иметь 5 по диагонали, а не 30, так как 30 был частью того же столбца, что и 10.

Я Я знаю, что для этого существует алгоритм под названием Венгерский алгоритм, но я не могу использовать для этой задачи никаких пакетов, кроме lpSolve.

Answer 1

Перестановка столбцов для матрицы A соответствует умножению матрицы AP, где P - матрица перестановок (переставленная единичная матрица). Таким образом, мы можем сформулировать следующую математическую модель:

Первое ограничение - Y=AP. Ограничения на P гарантируют, что P является правильной матрицей перестановок (по одному в каждой строке и столбце). Цель максимизирует след перестановочной матрицы столбца Y (след матрицы представляет собой сумму ее диагональных элементов).

Обратите внимание, что мы можем оптимизировать эту формулировку совсем немного (все y[i,j] с i<>j не используются, и мы можем заменить остальные у).

Некоторый код R, чтобы попробовать это:

library(CVXR)

# random matrix A
set.seed(123)
n <- 10
A <- matrix(runif(n^2,min=-1,max=1),nrow=n,ncol=n)

# decision variables
P <- Variable(n,n,boolean=T)
Y <- Variable(n,n)

# optimization model
# direct translation of the mathematical model given above
problem <- Problem(Maximize(matrix_trace(Y)),
                   list(Y==A %*% P,
                        sum_entries(P,axis=1) == 1,
                        sum_entries(P,axis=2) == 1))

# solve and print results
result <- solve(problem)
cat("status:",result$status)
cat("objective:",result$value)

В этом примере мы начинаем с матрицы

             [,1]        [,2]        [,3]        [,4]       [,5]       [,6]       [,7]        [,8]       [,9]       [,10]
 [1,] -0.42484496  0.91366669  0.77907863  0.92604847 -0.7144000 -0.9083377  0.3302304  0.50895032 -0.5127611 -0.73860862
 [2,]  0.57661027 -0.09333169  0.38560681  0.80459809 -0.1709073 -0.1155999 -0.8103187  0.25844226  0.3361112  0.30620385
 [3,] -0.18204616  0.35514127  0.28101363  0.38141056 -0.1725513  0.5978497 -0.2320607  0.42036480 -0.1647064 -0.31296706
 [4,]  0.76603481  0.14526680  0.98853955  0.59093484 -0.2623091 -0.7562015 -0.4512327 -0.99875045  0.5763917  0.31351626
 [5,]  0.88093457 -0.79415063  0.31141160 -0.95077263 -0.6951105  0.1218960  0.6292801 -0.04936685 -0.7942707 -0.35925352
 [6,] -0.90888700  0.79964994  0.41706094 -0.04440806 -0.7223879 -0.5869372 -0.1029673 -0.55976223 -0.1302145 -0.62461776
 [7,]  0.05621098 -0.50782453  0.08813205  0.51691908 -0.5339318 -0.7449367  0.6201287 -0.24036692  0.9699140  0.56458860
 [8,]  0.78483809 -0.91588093  0.18828404 -0.56718413 -0.0680751  0.5066157  0.6247790  0.22554201  0.7861022 -0.81281003
 [9,]  0.10287003 -0.34415856 -0.42168053 -0.36363798 -0.4680547  0.7900907  0.5886846 -0.29640418  0.7729381 -0.06644192
[10,] -0.08677053  0.90900730 -0.70577271 -0.53674843  0.7156554 -0.2510744 -0.1203366 -0.77772915 -0.6498947  0.02301092

Это имеет trace(A)=0.7133438.

Переменные Y имеют переставленные столбцы:

             [,1]        [,2]        [,3]        [,4]        [,5]        [,6]       [,7]       [,8]       [,9]      [,10]
 [1,]  0.92604847 -0.73860862  0.50895032  0.77907863 -0.42484496  0.91366669 -0.5127611  0.3302304 -0.9083377 -0.7144000
 [2,]  0.80459809  0.30620385  0.25844226  0.38560681  0.57661027 -0.09333169  0.3361112 -0.8103187 -0.1155999 -0.1709073
 [3,]  0.38141056 -0.31296706  0.42036480  0.28101363 -0.18204616  0.35514127 -0.1647064 -0.2320607  0.5978497 -0.1725513
 [4,]  0.59093484  0.31351626 -0.99875045  0.98853955  0.76603481  0.14526680  0.5763917 -0.4512327 -0.7562015 -0.2623091
 [5,] -0.95077263 -0.35925352 -0.04936685  0.31141160  0.88093457 -0.79415063 -0.7942707  0.6292801  0.1218960 -0.6951105
 [6,] -0.04440806 -0.62461776 -0.55976223  0.41706094 -0.90888700  0.79964994 -0.1302145 -0.1029673 -0.5869372 -0.7223879
 [7,]  0.51691908  0.56458860 -0.24036692  0.08813205  0.05621098 -0.50782453  0.9699140  0.6201287 -0.7449367 -0.5339318
 [8,] -0.56718413 -0.81281003  0.22554201  0.18828404  0.78483809 -0.91588093  0.7861022  0.6247790  0.5066157 -0.0680751
 [9,] -0.36363798 -0.06644192 -0.29640418 -0.42168053  0.10287003 -0.34415856  0.7729381  0.5886846  0.7900907 -0.4680547
[10,] -0.53674843  0.02301092 -0.77772915 -0.70577271 -0.08677053  0.90900730 -0.6498947 -0.1203366 -0.2510744  0.7156554

У нас есть trace(Y)=7.42218. Это лучшее, что мы можем сделать (доказано).

Answer 2

Это метод грубой силы, рассматривающий все перестановки. Это может стать несостоятельным для больших матриц.

library(RcppAlgos)
n = 5L
set.seed(123L)

mat = matrix(sample(1:10, n^2, TRUE), ncol = n)
mat
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    3    5    5    3    9
#> [2,]    3    4    3    8    3
#> [3,]   10    6    9   10    4
#> [4,]    2    9    9    7    1
#> [5,]    6   10    9   10    7

col_perms = permuteGeneral(n, n)
rows = seq_len(n)

diag_sum = apply(col_perms, 1, function(col) sum(mat[cbind(rows, col)]))
optim_cols = which.max(diag_sum)

mat[cbind(rows, col_perms[optim_cols, ])]
#> [1]  9  8 10  9 10
mat[, col_perms[optim_cols, ]]
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    9    3    3    5    5
#> [2,]    3    8    3    3    4
#> [3,]    4   10   10    9    6
#> [4,]    1    7    2    9    9
#> [5,]    7   10    6    9   10