Solve Transportation Problem by Aurenhammer–Hoffmann–Aronov Method


Solve transportation problem by Aurenhammer–Hoffmann–Aronov Method.


aha(a, b, nscales = 1, scmult = 2, factr = 1e+05, maxit = 10000, powerdiag=FALSE,
    wasser = FALSE, wasser.spt = NA, approx=FALSE, ...)
transport_apply(a, tplan)
transport_error(a, b, tplan)



an m x n matrix. a is treated as a measure on [0,m] x [0,n] with constant density on each unit square [i,i+1)x[j,j+1).


either a matrix such that dim(a)==dim(b) and sum(a)==sum(b) or a list of three vectors of equal length, named x, y and mass such that sum(a)==sum(b$mass), representing a discrete measure on [0,m]x[0,n].


a transference plan from a (to b), typically an optimal transference plan obtained by a call to aha.

nscales, scmult

the number of scales to use for the multiscale approach (the default is 1 meaning no multiscale approach), and the factor by which the number of pixels in each dimension is multiplied to get from a coarser to the next finer scale.

factr, maxit

parameters passed to the underlying L-BFGS-B algorithm (via the argument control in the R-function optim).


logical. Instead of an optimal transference plan, should the parameters for the optimal power diagram be returned?


logical. Instead of an optimal transference plan, should the L2-Wasserstein-distance between a and b be returned directly?


the number of support points used to approximate the discrete measure b. Defaults to NA meaning the full set of support points of b is used. If this argument is not NA, wasser is set to TRUE.


logical. If TRUE, an approximation to the objective function is used during optimization.


further arguments passed to optim via its argument control.


The function aha implements the algorithm by Aurenhammer, Hoffmann and Aronov (1998) for finding optimal transference plans in terms of the squared Euclidean distance in two dimensions. It follows the more detailed description given in Mérigot (2011) and also implements the multiscale version presented in the latter paper.

The functions transport_apply and transport_error serve for checking the accuracy of the transference plan obtained by aha. Since this transference plan is obtained by continuous optimization it will not transport exactly to the measure b, but to the measure transport_apply(a, tplan). By transport_error(a, b, tplan) the sum of absolut errors between the transported a-measure and the b-measure is obtained.


If powerdiag and wasser are both FALSE, a data frame with columns from, to and mass, which specify 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 and b. There are plot methods for the classes pgrid and pp, which can plot this solution.

If powerdiag is TRUE and wasser is FALSE, a list with components xi, eta, w and rect, which specify the parameters for the optimal power diagram in the same format as needed for the function power_diagram. Note that rect is always c(0,m,0,n).

If wasser is TRUE, a data frame with columns wasser.dist and error.bound of length one, where error.bound gives a bound on the absolute error in the Wasserstein distance due to approximating the measure b by a measure on a smaller number of support points.


Björn Bähre
(slightly modified by Dominic Schuhmacher


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

Q. Mérigot (2011). A multiscale approach to optimal transport. Eurographics Symposium on Geometry Processing 30(5), 1583–1592.

See Also

transport, which is a convenient wrapper function for various optimal transportation algorithms.


res <- aha(random32a$mass, random32b$mass)
plot(random32a, random32b, res, lwd=0.75)

aha(random64a$mass, random64b$mass, nscales=3, scmult=5, wasser.spt=3000, approx=TRUE)