Maximise diagonal of matrix by permuting columns in R

3
votes

I'm working with a square matrix in R, we can call it mat, and would like to permute the columns (i.e. change their order) so as to maximise the sum of the diagonal elements. I want to do this via linear programming methods, i.e. relying on the optimization package lpSolve. Code solutions are of course appreciated but failing that, any help formulating it as a linear programming problem would be appreciated.

My question is similar to this one: Permute columns of a square 2-way contingency table (matrix) to maximize its diagonal. However, in that question, and others I have found on SO, it's considered sufficient to go row-wise maximising the diagonal element in that row. The problem is that something like

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

isn't going to work for me: you could have situations where a row has multiple equal maxima, or where (say) in row X you pick 11 on the diagonal rather than 10 but consequently in row X+1 you are forced to have 5 on the diagonal rather than 30, since 30 was part of the same column as the 10.

I'm aware there is an algorithm called the Hungarian Algorithm for doing this, but I can't use any packages for this challenge except lpSolve.

2
rmatrixoptimizationlinear-programminglpsolve

2 Answers

6
votes

A column permutation for the matrix A corresponds to a matrix-multiplication AP where P is a permutation matrix (a permuted identity matrix). So we can formulate the following mathematical model:

enter image description here

The first constraint is Y=AP. The constraints on P make sure P is a proper permutation matrix (one 1 in each row and column). The objective maximizes the trace of the column-permuted matrix Y (the trace of a matrix is the sum of its diagonal elements).

Note that we can optimize this formulation quite a bit (all y[i,j] with i<>j are not used and we can substitute out the remaining y's).

Some R code to try this out:

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)

In this example, we start with the matrix

             [,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

This has trace(A)=0.7133438.

The Y variables have the columns permuted:

             [,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

We have trace(Y)=7.42218. This is the best we can do (proven).

2
votes

This is brute force method looking at all of the permutations. It's likely to become untenable for large matrices.

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