Ravi Varadhan <rvaradhan <at> jhmi.edu> writes:
>
> 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:
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
>
> 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
>
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
