I would like to minimize || PG-I || over P, where P is a p x p permutation matrix

(obtained by permuting the rows and/or columns of the identity matrix), G is a

given p x p matrix with full rank and I the identity matrix.

||.|| is the frobenius norm.

This can be implemented as a MIQP:

seti /i1*i25/;alias(i,j,k);

parameterG(i,j);

G(i,j) = uniform(-1,1);

parameterId(i,j);

Id(i,i) = 1;

binaryvariablep(i,j);variablez,v(i,j);

equation

defv(i,j)

norm

assign1(i)

assign2(j)

;

* Frobenius norm: http://mathworld.wolfram.com/FrobeniusNorm.html* we can drop the sqrt as this is a monotonic increasing function

norm.. z =e=sum((i,j), sqr(v(i,j)));

* we can substitute out v if needed

defv(i,j).. v(i,j) =e=sum(k, p(i,k)*g(k,j)) - Id(i,j);

* these assignment constraints make sure P is a permutation matrix

assign1(i)..sum(j, p(i,j)) =e= 1;

assign2(j)..sum(i, p(i,j)) =e= 1;

modelm /all/;solvem minimizing z using miqcp;

This is not very easy to solve to optimality if set i is larger than 20. The bounds just move too slowly. When running the problem for 1000 seconds with Cplex I see:

Although new solutions are found during the search even close to the end of the run, they only show marginal improvements.

Someone suggested to solve this as an assignment problem. I don’t believe that is actually 100% correct, but the approximation turns out to be very good. The GAMS code is:

parameterc(i,j);

c(i,j) = sqrt(sum(k, sqr(Id(i,k)-G(j,k))));

equationcassign;

cassign.. z =e=sum((i,j),c(i,j)*p(i,j));

modelassign/assign1,assign2,cassign/;solveassign minimizing z using rmip;

Here is small example that shows some differences:

seti /i1*i6/;alias(i,j,k);

parameterG(i,j);*G(i,j) = uniform(0,2);

G(i,j) = normal(0,2);

parameterId(i,j);

Id(i,i) = 1;

binaryvariablep(i,j);variablez,v(i,j);

equation

defv(i,j)

norm

assign1(i)

assign2(j)

;

* Frobenius norm: http://mathworld.wolfram.com/FrobeniusNorm.html* we can drop the sqrt as this is a monotonic increasing function

norm.. z =e=sum((i,j), sqr(v(i,j)));

* we can substitute out v if needed

defv(i,j).. v(i,j) =e=sum(k, p(i,k)*g(k,j)) - Id(i,j);

* these assignment constraints make sure P is a permutation matrix

assign1(i)..sum(j, p(i,j)) =e= 1;

assign2(j)..sum(i, p(i,j)) =e= 1;

modelm /all/;optionoptcr=0;solvem minimizing z using miqcp;

scalarfrobnorm1; frobnorm1 = sqrt(z.l);optionp:0;displayp.l,frobnorm1;

parameterc(i,j);

c(i,j) = sqrt(sum(k, sqr(Id(i,k)-G(j,k))));

equationcassign;

cassign.. z =e=sum((i,j),c(i,j)*p(i,j));

modelassign/assign1,assign2,cassign/;solveassign minimizing z using rmip;

scalarfrobnorm2,diff;

frobnorm2=sqrt(sum((i,j),sqr(sum(k, p.l(i,k)*g(k,j)) - Id(i,j))));

diff = frobnorm2-frobnorm1;displayp.l,frobnorm2,frobnorm1,diff;

The differences are:

---- 81 PARAMETER frobnorm2 = 10.584 |

A different suggestion seems to actually work. When we use

c(i,j) =

sum(k, sqr(Id(i,k)-G(j,k)));

the results with some random G matrices are the same. So to understand this we need to show that

give the same solutions.