通过对R列进行换列，使矩阵的对角线最大化。

Question

我在R中处理的是一个方形矩阵，我们可以把它叫做 mat我想对这些列进行换列（即改变它们的顺序），以使对角线元素的总和最大化。我想通过线性编程方法来实现这一目标，即依靠优化包lpSolve来实现。当然，如果能提供代码解决方案，我们将不胜感激，但如果不能，我们将感谢任何帮助，把它作为一个线性编程问题。

我的问题和这个问题类似。将一个正方形的双向应急表（矩阵）的列进行分解，使其对角线最大化。. 然而，在这个问题中，以及我在SO上发现的其他问题中，人们认为只要将该行的对角线元素最大化就可以了。问题是，像

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 是一个适当的换元矩阵（每行每列一个1）。目标最大限度地提高列permuted矩阵的轨迹 Y 矩阵的轨迹是其对角线元素的和）。

请注意，我们可以对这个公式进行相当程度的优化（所有的 y[i,j] 与 i<>j 不用，我们可以用剩下的y来代替）。)

一些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