Skip to content

optimization problem

16 messages · Erwin Kalvelagen, Ravi Varadhan, Klaus Nordhausen +1 more

Klaus Nordhausen
# 
Dear R-experts,

this is not a direct R-problem but I hope you can help me anyway.

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.

Does anyone know an algorithm to solve such a problem? And if yes, is it implemented in R?

Currently I minimize it by going through all possible permutations - but this is not feasible for higher dimensional problems as I would like to consider too.

Any help is appreciated!

Thanks in advance,

Klaus
Erwin Kalvelagen
# 
<klausch <at> gmx.de> writes:
matrix (obtained by permuting the rows
rank and I the identity matrix.
implemented in R?
is not feasible for higher
This could be modeled as a MIQP problem:

min sum((i,j),sqr(V(i,j)))
v(i,j) = sum(k, P(i,k)*G(k,j)) - Id(i,j);
sum(i, P(i,j)) = 1
sum(j, P(i,j)) = 1
P(i,j) in {0,1}

If you have access to Cplex then http://cran.r-
project.org/web/packages/Rcplex/Rcplex.pdf can help.

If you can use a different norm it may be possible to use linear MIP technology 
allowing a wider range of solvers.

This is combinatorial, so for p > 20 say it may become slow (this also depends 
on the data).

Erwin

----------------------------------------------------------------
Erwin Kalvelagen
Amsterdam Optimization Modeling Group
erwin at amsterdamoptimization.com
http://amsterdamoptimization.com
Ravi Varadhan
# 
Hi Klaus,

This problem can be cast as a linear sum assignment problem (LSAP), and solved using the `solve_LSAP' function in the "clue" package.  I show how to solve a slightly more general problem of finding a permulation of matrix A (n x n) that minimizes the Frobenius distance to a given matrix B (n x n).  In your case, B = I.  I ran this for n up to 100.  It seems to be quite fast.  

Here is how you do this:

dist <- function(A, B) { 
# Frobenius norm of A - B	
	n <- nrow(A)
	sum(abs(B - A))
}

pMatrix.min <- function(A, B) {
# finds the permutation P of A such that ||PA - B|| is minimum
# in Frobenius norm
# Uses the linear-sum assignment problem (LSAP) solver
# in the "clue" package
# Returns P%*%A and the permutation vector `pvec' such that
# A[pvec, ] is the permutation of A closest to B
	n <- nrow(A)
	D <- matrix(NA, n, n)
	for (i in 1:n) {
	for (j in 1:n) {
	D[j, i] <- sum(abs(B[j, ] - A[i, ]))
	} }
vec <- c(solve_LSAP(D))
list(A=A[vec,], pvec=vec)
}

#set.seed(123)
# Find P such that ||PA - I|| is minimum
n <- 100
A <- matrix(rnorm(n*n), n, n)
dist(A, diag(n))
B <- pMatrix.min(A, diag(n))$A  # B = P%*%A
dist(B, diag(n))

# Find P such that ||PA - C|| is minimum
C <- matrix(rnorm(n*n), n, n)
B <- pMatrix.min(A, C)$A
dist(A, C)
dist(B, C)

Let me know how this works for you.

Best,
Ravi.  

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: klausch at gmx.de
Date: Friday, January 15, 2010 9:53 am
Subject: [R] optimization problem
To: r-help at r-project.org
Erwin Kalvelagen
# 
Ravi Varadhan <rvaradhan <at> jhmi.edu> writes:
See http://mathworld.wolfram.com/FrobeniusNorm.html for a definition of the 
Frobenius norm.


Erwin

----------------------------------------------------------------
Erwin Kalvelagen
Amsterdam Optimization Modeling Group
erwin at amsterdamoptimization.com
http://amsterdamoptimization.com
Ravi Varadhan
# 
Thanks, Erwin, for pointing out this mistake.

Here is the correct function for Frobenius norm.  

Klaus - Just replace the old `dist' with the following one. 

dist <- function(A, B) {
 # Frobenius norm of A - B
  n <- nrow(A)
  sqrt(sum((B - A)^2))
 }

Ravi.

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Erwin Kalvelagen <erwin.kalvelagen at gmail.com>
Date: Saturday, January 16, 2010 2:35 am
Subject: Re: [R] optimization problem
To: r-help at stat.math.ethz.ch
Ravi Varadhan
# 
Klaus,

You also need to make a change in the main function, as shown below.

pMatrix.min <- function(A, B) {
# finds the permutation P of A such that ||PA - B|| is minimum
# in Frobenius norm
# Uses the linear-sum assignment problem (LSAP) solver
# in the "clue" package
# Returns P%*%A and the permutation vector `pvec' such that
# A[pvec, ] is the permutation of A closest to B
	n <- nrow(A)
	D <- matrix(NA, n, n)
	for (i in 1:n) {
	for (j in 1:n) {
#	D[j, i] <- sum(abs(B[j, ] - A[i, ]))
	D[j, i] <- sqrt(sum((B[j, ] - A[i, ])^2)) # correct Frobenius norm
	} }
vec <- c(solve_LSAP(D))
list(A=A[vec,], pvec=vec)
}

Hope this help,
Ravi.
____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Ravi Varadhan <rvaradhan at jhmi.edu>
Date: Saturday, January 16, 2010 10:00 am
Subject: Re: [R] optimization problem
To: Erwin Kalvelagen <erwin.kalvelagen at gmail.com>
Cc: r-help at stat.math.ethz.ch
Ravi Varadhan
# 
Yes, it can be formulated as an LSAP.  It works just fine.  I have checked it with several 3 x 3 examples.

Here is another convincing example:

n <- 50

A <- matrix(rnorm(n*n), n, n)
vec <- sample(1:n, n, rep=FALSE)  # a random permutation
  
C <- A[vec, ]  # the target matrix is just a permutation of original matrix
 
B <- pMatrix.min(A, C)$A
[1] 0
[1] 69.60859

Ravi.
____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Erwin Kalvelagen <erwin.kalvelagen at gmail.com>
Date: Saturday, January 16, 2010 1:36 pm
Subject: Re: [R] optimization problem
To: Ravi Varadhan <rvaradhan at jhmi.edu>
Cc: r-help at stat.math.ethz.ch
Ravi Varadhan
# 
Interesting! 

Now, if I change the "cost matrix", D,  in the LSAP formulation slightly such that it is quadratic, it finds the best solution to your example:


pMatrix.min <- function(A, B) {
# finds the permutation P of A such that ||PA - B|| is minimum
# in Frobenius norm
# Uses the linear-sum assignment problem (LSAP) solver
# in the "clue" package
# Returns P%*%A and the permutation vector `pvec' such that
# A[pvec, ] is the permutation of A closest to B
	n <- nrow(A)
	D <- matrix(NA, n, n)
	for (i in 1:n) {
	for (j in 1:n) {
#	D[j, i] <- sqrt(sum((B[j, ] - A[i, ])^2))
	D[j, i] <- (sum((B[j, ] - A[i, ])^2))  # this is better
	} }
vec <- c(solve_LSAP(D))
list(A=A[vec,], pvec=vec)
}

 > X<-pMatrix.min(A,B)
[1] 6 1 3 2 4 5
[1] 10.50172
This should be fine.  Any counter-examples to this?!

Best,
Ravi.
____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Erwin Kalvelagen <erwin.kalvelagen at gmail.com>
Date: Saturday, January 16, 2010 5:26 pm
Subject: Re: [R] optimization problem
To: Ravi Varadhan <rvaradhan at jhmi.edu>
Cc: r-help at stat.math.ethz.ch
Hans W Borchers
# 
Ravi Varadhan <rvaradhan <at> jhmi.edu> writes:
Dear Ravi,

I thought your solution is ingenious, but after the discussion with 
Erwin Kalvelagen I found two things quite irritating:

(1) Why is solve_LSAP(D) different from solve_LSAP(D^2) in Erwin's
    example? I believed just squaring the distance matrix would make
    no difference to solving the LSAP - except for some numerical
    instability which does not seem to be the case here.

(2) If you change rows and sums in the definition of D, that is

    D[j, i] <- sqrt(sum((B[, j] - A[, i])^2))

    then the solution to Erwin's example comes out right even with
    keeping the square root.

Do you have explanations for these 'phenomena'? Otherwise, I think,
there will remain some doubts about this approach.

Very best
Hans Werner
Klaus Nordhausen
# 
Dear Erwin, Ravi and Hans Werner,

thanks a lot for your replies. I don't think I have access to Cplex and therefore probably cannot try that out but will read about MIQP.

I played a bit then around with Ravi's suggestions and made also the observation that the linear cost function often found the exact solution but now always - but in my tests the quadratic cost function version always found the correct result. But I'll test that still a bit more detailed.

Would anyone of you know a good reference for an overview what algorithms are there and for which problems they can be used?

Thank a lot again!

Klaus

-------- Original-Nachricht --------

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ravi Varadhan
    

    

      
      #
    

  


  

      
Dear Hans,

I agree with your comments.  My intuition was that the quadratic form would be better behaved  than the radical form (less nonlinear!?).  So, I was "hoping" to see a change in behavior when the cost function was altered from a radical (i.e. sqrt) form to quadratic, but I was still surprised to actually see it.  I don't have a good explanation for this.  

I came up with the idea of solving Klaus' problem as an LSAP problem.  I did not know that the LSAP approach to this and similar problems has already been considered by Nick Higham.  I asked Nick last night about this problem thinking that he might know of more direct solutions to this problem (e.g. some kind of SVD or related factorizations). He said that I should take a look at the PhD thesis of one of his students.

Take a look at Section 3.5 of the PhD thesis

   Parallel Solution of SVD-Related Problems, With Applications

at

http://www.maths.manchester.ac.uk/~higham/misc/past-students.php


This thesis proposed algorithms for the current problem and different versions of it under the heading of "Procrustes-type" problems.  I have to read this and get a better handle on it.  However, I would not be able to get to this for another two weeks.  If you have any insights from this, in the meanwhile, do share with us.

Best regards,
Ravi.

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: "Hans W. Borchers" <hwborchers at googlemail.com>
Date: Sunday, January 17, 2010 3:54 am
Subject: Re: [R] optimization problem
To: r-help at stat.math.ethz.ch

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Ravi Varadhan
    

    

      
      #
    

  


  

      
Klaus,

I am happy to know that the quadratic cost LSAP seems to work well for you.  

The Hungarian algorithm is a classic for solving linear sum assignment problem, which is closely related to matching in bipartite graphs.  You can google or wiki these terms to get papers and books on this topic. 

Also look at the PhD thesis by one of Nick Higham's students (Section 3.5):

 Parallel Solution of SVD-Related Problems, With Applications

at

http://www.maths.manchester.ac.uk/~higham/misc/past-students.php

This thesis has a good discussion of different variants of your problem, some of which are more general than yours.  

Your problem is one of the variants of "procrustes-type" problems found in multivariate statistics.  The very same problem that you posed is supposed to occur in multidimensional scaling.  So, you might also want to look in that literature.


Best,
Ravi.

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: klausch at gmx.de
Date: Sunday, January 17, 2010 8:06 am
Subject: Re: [R] optimization problem
To: Ravi Varadhan <rvaradhan at jhmi.edu>, erwin.kalvelagen at gmail.com, hwborchers at googlemail.com
Cc: r-help at stat.math.ethz.ch

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Hans W Borchers
    

    

      
      #
    

  


  

      
Ravi Varadhan <rvaradhan <at> jhmi.edu> writes:

            

      
      
      
    

        
Thanks for pointing out the thesis. As I understand, the (one-sided)
Procrustes problem is finding an orthogonal matrix minimizing

    min! || A R - B || , t(R) R = I , ||.|| the Frobenius norm

and that there is an substantial theory behind in numerical linear
algebra. The basic literature appears to be

    Gower, J.C; Dijksterhuis, G.B., Procrustes Problems, Oxford
    Statistical Science Series, Oxford University Press, 2004

The thesis itself will lead us astray as it treats the two-sided
Procrustes Problem, which is not our main concern.
The request on R-help only asked for permutation matrices P (from
left or right) minimizing

    min! || P A - B ||

so I still think that a direct approach as you have suggested is
possible given this apparently much simpler problem.

Take your definition of D with quadratic terms:
The LSAP approach will find a permutations of the rows of B, say Bp,
such that the (linear!) sum over D_{ip(i)} is minimal, that is

    sum over rows {sum d(Bp - A)^2} = sum over all {(Bp - A)^2}

which is exactly square of the Frobenius norm.
If you instead apply your first definition with square roots, it is

    sum over rows {sum sqrt(d(Bp - A)^2)}

and this cannot be expanded to the sum of the Frobenius norm.
Therefore, the quadratic approach is indeed correct and will lead
to a true optimum, while the first variant will not.

Hans Werner