Integer nonlinear optimization with a fixed domain

Question

Here is the original table:

       Column1   Column2
Row1    4 (x1)    6(x3)
Row2    5 (x2)    4(x4)

x1, x2, x3, x4 mean the new values after swapping

Originally the variance of each column is:

Var(column1)= var(c(4,5))=0.5
Var(column2)= var(c(6,4))=2

After swapping the values lesser than 10 in the original table, the new variances of two columns respectively are:

New_Var(column1)=(x1-(x1+x2)/2)^2+(x2-(x1+x2)/2)^2
New_Var(column2)=(x3-(x3+x4)/2)^2+(x4-(x3+x4)/2)^2

My objective is to

Minimize  | New_Var(column1) – 0.5 | + | New_Var(column2) – 2 |

Here the notation ‘| |’ means absolute.

The constrains are that x1, x2, x3, x4 can only get a value respectively from a fixed domain, here say {4, 5, 4, 6}. And all the four variables need be swapped to a different value, say that x1 cannot be 4, x3 cannot be 6 so on so forth.

Note that all values in the original table here are larger than 10 because I just want to make the question look like simple.

In realistic, the table is 6000*10 which is very large. So outputting all the unique permutations and testing is not a suitable method.

I have read the task view of optimization in R. There are a lot of optimization packages there. So I need more specific guidance.

Below is the function that can swap values lesser than 10 to a different values lesser than 10 respectively. Hope it helps.

derangement <- function(x){
  if(max(table(x)) > length(x)/2) return(NA)
  while(TRUE){
    y <- sample(x)
    if(all(y != x)) return(y)
  }
}

swapFun <- function(x, n = 10){
  inx <- which(x < n)
  y <- derangement(x[inx])
  if(length(y) == 1) return(NA) 
  x[inx] <- y
  x
}

set.seed(10)
swapFun(c(1,2,3,10,4,11,2,12))
#[1]  2  4 10  2 11  1 12

How can I solve it? Thanks.

Answer 1

For problems this size you can easily use brute force.

In R:

library(combinat)
pp <- permn(c(1,2,2,3)) ## construct list of permutations
pp <- unique(pp)        ## unique permutations only

Apply your (new) constraint "x1 cannot be 1, x2 cannot be 2, x3 cannot be 2, x4 cannot be 3":

ok <- function(X) {
   !(X[1]==1 | X[2]==2 | X[3]==2 | X[4]==3)
}
pp <- pp[sapply(pp,ok)]

There are length(pp) == 2 non-identical, valid combinations (you can test whether length(pp)>0 to see if there are any valid combinations left).

## define objective function
f <- function(X) { (X[1]-X[2])^2 + (X[1]-X[3])^2+ 
        (X[2]-X[3]-X[4])^2+X[4]^2+X[3]^2 }
## compute obj function for all permutations
vals <- sapply(pp,f)
## identify minima
pp[which(vals==min(vals))]
## [[1]]
## [1] 2 3 1 2 

This will break down as the domain grows and the size of (n!) gets large ... at which point you may need to start doing research in the domain of integer programming (e.g. see the Optimization Task View)

In Python you would first generate the list of permutations with itertools.permutations (see this question) and then use a comprehension to apply the objective function across the iterator.

Answer 2

I write this sometimes as a set of linear assignment constraints:

 sum(i, δ[i,j]) = 1  forall j
 sum(j, δ[i,j]) = 1  forall i
 x[i] = sum(j, δ[i,j] * v[j])
 δ[i,j] ∈ {0,1}     (binary variables)
 v = [1, 2, 2, 3]   (constants)  
 i,j ∈ {1..4}       (indices} 

Here δ is a binary variable that performs like a permutation matrix.