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 key word used. Much more details can be found in the cited references
and in a forthcoming package vignette.
aha: The Aurenhammer--Hoffmann--Aronov (1998) method with the multiscale approach presented in Mérigot (2011). The original theory was limited to \(p=2\). We refer by aha also to the extension of the same idea for \(p=1\) as presented in Hartmann and Schuhmacher (2017) and for more general \(p\) (currently not implemented).
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).
networkflow: The fast implementation of the network simplex algorithm by Nicolas Bonneel based on the LEMON Library (see citations below).
primaldual: The primal-dual algorithm as described in Luenberger (2003, Section 5.9).
revsimplex: The revised simplex algorithm as described in Luenberger and Ye (2008, Section 6.4) with various speed improvements, including a multiscale approach.
shielding: The shielding (or shortcut) method, as described in Schmitzer (2016).
shortsimplex: The shortlist method based an a revised simplex algorithm, as described in Gottschlich and Schuhmacher (2014).
The order of the default key words specified for the argument method gives a rough idea of the relative efficiency of the algorithms for the corresponding class of objects. For a given a and b the actual computation times may deviate significantly from this order.
For class pgrid the default method is "auto", which resolves to "revsimplex" if p is not 2 or the problem is very small, and to "shielding" otherwise.
The following table gives information about the applicability of the various algorithms (or sometimes rather
their current implementations).
| default | pgrid | pp | wpp |
| aha (p=1 or 2!) | - | + | - | @ |
| auction | - | - | + | - |
| auctionbf | - | - | + | - |
| networkflow | + | + | + | + |
| primaldual | * | * | * | + |
| revsimplex | + | + | * | + |
| shielding (p=2!) | - | + | - | - |
| shortsimplex | + | + | * | + |
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 |
aha (\(p=2\)!) | - | + | factr, maxit |
auction | - | - | lasteps, epsfac |
auctionbf | - | - | lasteps, epsfac |
networkflow | - | - | |
primaldual | - | - | |
revsimplex | + | + | |
shielding (\(p=2\)!) | - | + | |
shortsimplex | - | - | slength, kfound, psearched |
start specifies the algorithm for computing a starting solution (if needed). 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. However,
for transport.pgrid and p larger than 1, there are two cases where an automatic multiscale procedure is also performed, i.e. the optimal transport is first computed on coarser grids and information from these solutions is then used for the finer girds.
This happens for
method = "revsimplex", where a single coarsening at factor scmult=2 is performed, and for method = "shielding", where a number of coarsenings adapted to the dimensions of the array is performed.
For p=1 and method="revsimplex", as well as p=2 and method="aha" there are multiscale versions of
the corresponding algorithms that allows for finer control via the parameters
nscales, scmult and returncoarse. 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).
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.
slength, kfound, psearched are the shortlist length, the number of pivot candidates needed, and the percentage of
shortlists searched, respectively.