# Find Optimal Transference Plan (a.k.a. Optimal Transport Map) Between Two Objects.

### Description

Given two objects `a`

and `b`

that specify distributions of mass and an object that specifies (a way to compute) costs,
find the transference plan for going from `a`

to `b`

that minimizes the total cost.

### Usage

1 2 3 4 5 6 7 8 9 10 11 12 13 | ```
transport(a, b, ...)
## Default S3 method:
transport(a, b, costm, method = c("shortsimplex", "revsimplex", "primaldual"),
control = list(), ...)
## S3 method for class 'pgrid'
transport(a, b, p = NULL, method = c("revsimplex", "shortsimplex",
"primaldual", "aha"), control = list(), ...)
## S3 method for class 'pp'
transport(a, b, p = 1, method = c("auction", "auctionbf", "shortsimplex",
"revsimplex", "primaldual"), control = list(), ...)
## S3 method for class 'wpp'
transport(a, b, p = 1, method = c("revsimplex", "shortsimplex", "primaldual"),
control = list(), ...)
``` |

### Arguments

`a, b` |
two objects that describe mass distributions, between which the optimal transport map is to be computed. For the default
method these are vectors of non-negative values. For the other three methods these are objects of the respective classes.
It is also possible to have |

`costm` |
for the default method a |

`p` |
for the three specialized methods the power |

`method` |
the name of the algorithm to use. See details below. |

`control` |
a named list of parameters for the chosen method or the result of a call to |

`...` |
currently without effect. |

### Details

There is a number of algorithms that are currently implemented and more will be added in future versions of the package. The following is a brief description of each abbreviation used. Much more details can be found in the cited references and in a forthcoming package vignette.

`shortsimplex`

: The shortlist method based an a revised simplex algorithm, as described in Gottschlich and Schuhmacher (2014).

`revsimplex`

: The revised simplex algorithm as described in Luenberger and Ye (2008, Section 6.4) with various speed improvements, including a multiscale approach.

`primaldual`

: The primal-dual algorithm as described in Luenberger (2003, Section 5.9).

`aha`

: The Aurenhammer–Hoffmann–Aronov (1998) method with the multiscale approach presented in Mérigot (2011).

`auction`

: The auction algorithm by Bertsekas (1988) with epsilon-scaling, see Bertsekas (1992).

`auctionbf`

: A refined auction algorithm that combines forward and revers auction, see Bertsekas (1992).

The following table gives information about the applicability of the various algorithms (or sometimes rather their current implementations).

default | pgrid | pp | wpp | |

shortsimplex | + | + | * | + |

revsimplex | + | + | * | + |

primaldual | * | * | * | + |

aha (p=2!) | - | + | - | @ |

auction | - | - | + | - |

auctionbf | - | - | + | - |

where: + recommended, * applicable (may be slow), - no implementation planned or combination does not make sense; @ indicates that the aha algorithm is available in the special combination where `a`

is a `pgrid`

object and `b`

is a `wpp`

object (and `p`

is 2). For more details on this combination see the function `semidiscrete`

.

Each algorithm has certain parameters supplied by the `control`

argument. The following table gives an overview of parameter names and
their applicability.

`start` | multiscale | individual parameters | |

`shortsimplex` | - | - | `slength` , `kfound` , `psearched` |

`revsimplex` | + | + | |

`primaldual` | - | - | |

`aha` (p=2!) | - | + | `factr` , `maxit` |

`auction` | - | - | `lasteps` , `epsfac` |

`auctionbf` | - | - | `lasteps` , `epsfac` |

`start`

specifies the algorithm for computing a starting solution. Currently the Modified Row Minimum Rule
(`start="modrowmin"`

), the North-West Corner Rule (`start="nwcorner"`

) and the method by Russell (1969) (`start="russell"`

)
are implemented. When `start="auto"`

(the default) the ModRowMin Rule is chosen, except in `transport.pgrid`

if `p`

>1,
where the optimal transport is first computed on a grid that is coarser by a factor of `scmult=2`

and a refinement of this solution
is used as a starting solution.

For `p=1`

and `method="revsimplex"`

, as well as `p=2`

and `method="aha"`

there is a proper multiscale versions of
the corresponding algorithm that allows for finer control over a similar behaviour and is based on separate code.
There are three parameters `nscales`

, `scmult`

, `returncoarse`

(summarized under the title
“multiscale” in the above table). The default value of `nscales=1`

suppresses
the multiscale version. For larger problems it is advisable to use the multiscale version, which currently is only implemented for
square pgrids in two dimensions. The algorithm proceeds then by coarsening the pgrid `nscales-1`

times, summarizing
each time `scmult^2`

pixels into one larger pixels, and then solving the various transport problems starting from the coarsest and
using each previous problem to compute a starting solution to the next finer problem. If `returncoarse`

is `TRUE`

, the coarser
problems and their solutions are returned as well (`revsimplex`

only).

`slength`

, `kfound`

, `psearched`

are the shortlist length, the number of pivot candidates needed, and the percentage of
shortlists searched, respectively.

`factr`

, `maxit`

are the corresponding components of the `control`

argument in the `optim`

L-BFGS-B method.

`lasteps`

, `epsfac`

are parameters used for epsilon scaling in the auction algortihm. The algorithm starts with a “transaction cost” per bid of `epsfac^k * lasteps`

for some reasonable `k`

generating finer and finer approximate solutions as the `k`

counts down to zero. Note that in order for the procedure to make sense, `epsfac`

should be larger than one (typically two- to three-digit) and in order for the final solution to be exact `lasteps`

should be smaller than `1/n`

, where `n`

is the total number of points in either of the point patterns.

### Value

A data frame with columns `from`

, `to`

and `mass`

that specifies from which knot to which other knot what amount of mass is sent in the optimal transference plan. Knots are given as indices in terms of the usual column major enumeration of the matrices `a$mass`

and `b$mass`

. There are `plot`

methods for the classes `pgrid`

and `pp`

, which can plot this solution.

If `returncoarse`

is `TRUE`

for the `revsimplex`

method, a list with components `sol`

and `prob`

giving the solutions and problems on the various scales considered. The solution on the finest scale (i.e. the output we obtain when setting `returncoarse`

to `FALSE`

) is in `sol[[1]]`

.

If `a`

is of class `pgrid`

and `b`

of class `wpp`

(and `p=2`

), an object of class `power_diagram`

as described in the help for the function `semidiscrete`

. The `plot`

method for class `pgrid`

can plot this solution.

### Author(s)

Dominic Schuhmacher dschuhm1@uni-goettingen.de

Björn Bähre bjobae@gmail.com

(code for `aha`

-method)

Carsten Gottschlich gottschlich@math.uni-goettingen.de

(original java code for `shortlist`

- and `revsimplex`

-methods)

### References

F. Aurenhammer, F. Hoffmann and B. Aronov (1998). Minkowski-type theorems and least-squares clustering. Algorithmica 20(1), 61–76.

D. P. Bertsekas (1988). The auction algorithm: a distributed relaxation method for the assignment problem. Annals of Operations Research 14(1), 105–123.

D. P. Bertsekas (1992). Auction algorithms for network flow problems: a tutorial introduction. Computational Optimization and Applications 1, 7–66.

C. Gottschlich and D. Schuhmacher (2014). The shortlist method for fast computation of the earth mover's distance and finding optimal solutions to transportation problems. PLOS ONE 9(10), e110214. doi:10.1371/journal.pone.0110214

D.G. Luenberger (2003). Linear and nonlinear programming, 2nd ed. Kluwer.

D.G. Luenberger and Y. Ye (2008). Linear and nonlinear programming, 3rd ed. Springer.

Q. Mérigot (2011). A multiscale approach to optimal transport. Computer Graphics Forum 30(5), 1583–1592. doi:10.1111/j.1467-8659.2011.02032.x

### See Also

`plot`

, `wasserstein`

### Examples

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 | ```
#
# example for the default method
#
a <- c(100, 200, 80, 150, 50, 140, 170, 30, 10, 70)
b <- c(60, 120, 150, 110, 40, 90, 160, 120, 70, 80)
set.seed(24)
costm <- matrix(sample(1:20, 100, replace=TRUE), 10, 10)
res <- transport(a,b,costm)
# pretty-print solution in matrix form for very small problems:
transp <- matrix(0,10,10)
transp[cbind(res$from,res$to)] <- res$mass
rownames(transp) <- paste(ifelse(nchar(a)==2," ",""),a,sep="")
colnames(transp) <- paste(ifelse(nchar(b)==2," ",""),b,sep="")
print(transp)
#
# example for class 'pgrid'
#
dev.new(width=9, height=4.5)
par(mfrow=c(1,2), mai=rep(0.1,4))
image(random32a$mass, col = grey(0:200/200), axes=FALSE)
image(random32b$mass, col = grey(0:200/200), axes=FALSE)
res <- transport(random32a,random32b)
dev.new()
par(mai=rep(0,4))
plot(random32a,random32b,res,lwd=1)
#
# example for class 'pp'
#
set.seed(27)
x <- pp(matrix(runif(1000),500,2))
y <- pp(matrix(runif(1000),500,2))
res <- transport(x,y)
dev.new()
par(mai=rep(0.02,4))
plot(x,y,res)
#
# example for class 'wpp'
#
set.seed(30)
m <- 30
n <- 60
massx <- rexp(m)
massx <- massx/sum(massx)
massy <- rexp(n)
massy <- massy/sum(massy)
x <- wpp(matrix(runif(2*m),m,2),massx)
y <- wpp(matrix(runif(2*n),n,2),massy)
res <- transport(x,y,method="revsimplex")
plot(x,y,res)
#
# example for semidiscrete transport between class 'pgrid' and class 'wpp'
#
set.seed(33)
n <- 100
massb <- rexp(n)
massb <- massb/sum(massb)*1e5
b <- wpp(matrix(runif(2*n),n,2),massb)
res <- transport(random32a,b,p=2)
plot(random32a,b,res)
``` |