[R] optimization problem

[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)
> X$pvec
[1] 6 1 3 2 4 5
> dist(X$A, B)
[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


> I believe this is a very good approximation but not a 100% correct
> formulation of the original problem.
> 
> E.g. for
> 
> 
> A <- matrix(c(
>    -0.62668585718,-0.78785288063,-1.32462887089, 0.63935994044,
>    -1.99878497801, 4.42667400292, 0.65534961645, 1.86914537669,
>    -0.97229929674, 2.37404268115, 0.01223810011,-1.24956493590,
>     0.92711756473,-1.51859351813, 3.76707054743,-2.30641777527,
>    -1.98918429361,-0.43634856903,-3.65989556798, 0.00073904070,
>    -1.44153837918, 1.42004180214,-1.81388322489,-1.92423923917,
>    -1.46322482727,-1.81783134835,-2.59801302663, 2.03462135096,
>     1.32546711550,-0.21519150247,-1.94319654347, 0.68773346973,
>    -2.75094791807,-1.44814080195, 3.14197123843,-0.52521369103
> ),nrow=6)
> 
> B<-diag(nrow=6)
> 
> X<-pMatrix.min(A,B)
> 
> dist(X$A, B)
> 
> X$pvec
> 
> bestpvec <- c(6,1,3,2,4,5)
> 
> dist(A[bestpvec,],B)
> 
> you get a norm of 10.58374 while the true optimal norm is 10.50172.  For
> small problems you often get the optimal solution, but the error 
> caused by
> linearizing the objective becomes larger if the problems are larger. 
> But the
> approximation is actually very good.
> 
> Erwin
> 
> 
> ----------------------------------------------------------------
> Erwin Kalvelagen
> Amsterdam Optimization Modeling Group
> erwin at amsterdamoptimization.com
> 
> ----------------------------------------------------------------
> 
> 
> On Sat, Jan 16, 2010 at 2:01 PM, Ravi Varadhan <rvaradhan at jhmi.edu> wrote:
> 
> >
> > 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)
> >
> > > # Find P such that ||PA - C|| is minimum
> >
> > 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
> >
> > > dist(B, C)
> > [1] 0
> >
> > > dist(A, C)
> > [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
> >
> >
> > > I also have doubts this can be formulated correctly as a linear
> > assignment
> > > problem. You may want to check the results with a small example.
> > >
> > > Erwin
> > >
> > > ----------------------------------------------------------------
> > > Erwin Kalvelagen
> > > Amsterdam Optimization Modeling Group
> > > erwin at amsterdamoptimization.com
> > >
> > > ----------------------------------------------------------------
> > >
> > >
> > > On Sat, Jan 16, 2010 at 9:59 AM, Ravi Varadhan <rvaradhan at jhmi.edu>
> > wrote:
> > >
> > > >
> > > > 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 <rvaradhan <at> jhmi.edu> writes:
> > > > > > dist <- function(A, B) {
> > > > > > # Frobenius norm of A - B
> > > > > >     n <- nrow(A)
> > > > > >     sum(abs(B - A))
> > > > > > }
> > > > > >
> > > > >
> > > > > See  for a definition of the
> > > > > Frobenius norm.
> > > > >
> > > > >
> > > > > Erwin
> > > > >
> > > > > ----------------------------------------------------------------
> > > > > Erwin Kalvelagen
> > > > > Amsterdam Optimization Modeling Group
> > > > > erwin at amsterdamoptimization.com
> > > > >
> > > > >
> > > > > ______________________________________________
> > > > > R-help at r-project.org mailing list
> > > > >
> > > > > PLEASE do read the posting guide
> > > > > and provide commented, minimal, self-contained, reproducible code.
> > > >
> >