Title: | Approximate Optimal Transport Between Two-Dimensional Grids |
---|---|
Description: | Can be used for optimal transport between two-dimensional grids with respect to separable cost functions of l^p form. It utilizes the Frank-Wolfe algorithm to approximate so-called pivot measures: one-dimensional transport plans that fully describe the full transport, see G. Auricchio (2021) <arXiv:2105.07278>. For these, it offers methods for visualization and to extract the corresponding transport plans and costs. Additionally, related functions for one-dimensional optimal transport are available. |
Authors: | Michel Groppe [aut, cre], Nicholas Bonneel [ctb], Egerváry Research Group on Combinatorial Optimization [cph] |
Maintainer: | Michel Groppe <[email protected]> |
License: | GPL-3 |
Version: | 1.0.1 |
Built: | 2024-11-12 06:34:37 UTC |
Source: | CRAN |
Calculate an optimal dual pair for the optimal transport between discrete one-dimensional measures.
dual1d(a, b, wa, wb, p = 1, right.margin = 1e-15, sorted = FALSE)
dual1d(a, b, wa, wb, p = 1, right.margin = 1e-15, sorted = FALSE)
a |
first vector of points. |
b |
second vector of points. |
wa |
weight vector of the first vector of points. |
wb |
weight vector of the second vector of points. |
p |
the power |
right.margin |
small amount the points are moved by. |
sorted |
logical value indicating whether or not |
The pair is an optimal dual pair if the optimal transport distance between the two distributions
with respect to the cost function
is given by
and the condition holds.
a list containing the dual vectors pot.a
and pot.b
.
set.seed(1) a <- 1:5 wa <- rep(1/5, 5) b <- 1:6 wb <- runif(6) wb <- wb / sum(wb) d <- dual1d(a, b, wa, wb, p = 1) dc <- sum(d$pot.a * wa) + sum(d$pot.b * wb) print(all.equal(dc, transport_cost(a, b, wa, wb, p = 1)))
set.seed(1) a <- 1:5 wa <- rep(1/5, 5) b <- 1:6 wb <- runif(6) wb <- wb / sum(wb) d <- dual1d(a, b, wa, wb, p = 1) dc <- sum(d$pot.a * wa) + sum(d$pot.b * wb) print(all.equal(dc, transport_cost(a, b, wa, wb, p = 1)))
Calculate the transport plan obtained by the north-west-corner rule.
north_west_corner(wx, wy, threshold = 1e-15)
north_west_corner(wx, wy, threshold = 1e-15)
wx |
first weight vector. |
wy |
second weight vector. |
threshold |
small value that indicates when a value is considered to be zero. |
a matrix representing the transport plan obtained by the north-west-corner rule.
set.seed(1) wx <- rep(1/5, 5) wy <- runif(6) wy <- wy / sum(wy) north_west_corner(wx, wy)
set.seed(1) wx <- rep(1/5, 5) wy <- runif(6) wy <- wy / sum(wy) north_west_corner(wx, wy)
Create an object that represents a probability measure that is supported on a two-dimensional grid.
otgrid( mass, x = NULL, y = NULL, sorted = FALSE, normalize = FALSE, remove.zeros = FALSE )
otgrid( mass, x = NULL, y = NULL, sorted = FALSE, normalize = FALSE, remove.zeros = FALSE )
mass |
matrix of non-negative weights. |
x |
vector of support points of the first marginal. |
y |
vector of support points of the second marginal. |
sorted |
logical value specifying whether or not the support points are sorted. |
normalize |
logical value specifying whether or not the total mass should be rescaled to 1. |
remove.zeros |
logical value specifying whether or not marginals with no mass should be removed from the grid. |
If x
and y
are not specified, then a equispaced unit grid is chosen.
object of class "otgrid"
. It contains the following elements:
x |
vector of support points of the first marginal |
y |
vector of support points of the second marginal |
n |
number of support points of the first marginal |
m |
number of support points of the second marginal |
mass |
matrix of non-negative weights |
total |
total mass |
Also note that the functions print
and plot
are available for objects of class "otgrid"
.
plot plot.otgrid
x <- otgrid(cbind(1:2, 0, 3:4), remove.zeros = TRUE) print(x) # note that it's only 2 x 2 plot(x)
x <- otgrid(cbind(1:2, 0, 3:4), remove.zeros = TRUE) print(x) # note that it's only 2 x 2 plot(x)
Calculate the pivot measure of the optimal transport between two-dimensional grids.
pivot_measure(x, ...) ## S3 method for class 'otgrid' pivot_measure( x, y, p.1 = 2, p.2 = p.1, max.it = 100, tol = 1e-04, threads = 1, start.pivot = c("independent", "northwestcorner"), dual.method = c("discrete", "epsilon-discrete", "epsilon-histogram"), dual.params = list(), return.it = FALSE, ... ) ## S3 method for class 'data.frame' pivot_measure(x, a, b, p.1 = 2, p.2 = p.1, ...)
pivot_measure(x, ...) ## S3 method for class 'otgrid' pivot_measure( x, y, p.1 = 2, p.2 = p.1, max.it = 100, tol = 1e-04, threads = 1, start.pivot = c("independent", "northwestcorner"), dual.method = c("discrete", "epsilon-discrete", "epsilon-histogram"), dual.params = list(), return.it = FALSE, ... ) ## S3 method for class 'data.frame' pivot_measure(x, a, b, p.1 = 2, p.2 = p.1, ...)
x |
an object of class |
... |
further arguments (currently unused). |
y |
an object of class |
p.1 |
the first power |
p.2 |
the second power |
max.it |
the maximum number of iterations. |
tol |
the desired accuracy. |
threads |
number of threads to use. |
start.pivot |
the start pivot to use: matrix representing an actual coupling or |
dual.method |
the name of the dual calculators to use: |
dual.params |
list of parameters for the dual calculators. |
return.it |
logical value specifying whether or not the costs and dual gaps in each iteration are to be returned . |
a |
an object of class |
b |
an object of class |
Denote with and
the two-dimensional grids that
x
and y
lie on, respectively.
The pivot measure is a measure on that specifies the whole transport between
x
and y
given that the cost
is of the separable form .
The pivot measure is approximated using the Frank-Wolfe algorithm. The algorithm starts with an initial guess (start.pivot
), e.g., the
independent coupling ("independent"
) or the north-west-corner rule ("northwestcorner"
). Then, in each iteration step, the dual solutions of
multiple one-dimensional transport problems are calculated and combined to give the gradient. To ensure differentiability, the dual solutions
must be unique. There are three different calculators for the dual solutions, of which the last two ensure uniqueness:
name | description | dual.params |
"discrete" |
the distributions are basically unchanged, due to rounding errors the support points are moved by a small positive amount. | right.margin = 1e-15 |
"epsilon-discrete" |
additionally, the support is made connected by uniformly adding
mass. |
+ eps = 1e-8 |
"epsilon-histogram" |
additionally, each point mass is distributed uniformly in a bin around said point | + width = 1e-8
|
A pivot measure can also be computed from an optimal transport plan.
Finally, the pivot measure can be used to calculate the full transport plan and cost.
an object of class "otgridtransport"
that contains the following elements:
from |
an object of class "otgrid" the mass is transported from. |
to |
an object of class "otgrid" the mass is transported to. |
p.1 |
the first power of the cost function. |
p.2 |
the second power of the cost function. |
pivot |
a matrix representing the pivot measure. |
If the Frank-Wolfe algorithm is used to approximate the pivot measure, the element conv
indicates whether or not
we can ensured that the error is less or equal to the given precision tol
.
If return.it = TRUE
, then it also contains the vectors costs
and dualgaps
of costs and dual gaps in each iteration.
Also note that the functions print
and plot
are available for objects of class "otgridtransport"
.
transport plan transport_df
, transport cost transport_cost
, two-dimensional grid otgrid
,
plot plot.otgridtransport
x <- otgrid(rbind(c(0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 1, 0, 0, 0, 0, 0, 0), c(0, 1, 2, 1, 0, 0, 0, 0, 0), c(0, 0, 1, 0, 0, 0, 0, 0, 0), 0, 0, 0, 0, 0)) y <- otgrid(rbind(0, 0, 0, 0, c(0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 0, 0, 0, 0, 1, 0, 0), c(0, 0, 0, 0, 0, 1, 2, 1, 0), c(0, 0, 0, 0, 0, 0, 1, 0, 0), 0)) pm <- pivot_measure(x, y, p.1 = 1, p.2 = 2, dual.method = "epsilon-discrete") print(pm) plot(pm) # use pivot measure to calculate cost and plan pm <- transport_cost(pm) pm <- transport_df(pm) print(pm) # calculate pivot from plan pm2 <- pivot_measure(pm$df, x, y) plot(pm2)
x <- otgrid(rbind(c(0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 1, 0, 0, 0, 0, 0, 0), c(0, 1, 2, 1, 0, 0, 0, 0, 0), c(0, 0, 1, 0, 0, 0, 0, 0, 0), 0, 0, 0, 0, 0)) y <- otgrid(rbind(0, 0, 0, 0, c(0, 0, 0, 0, 0, 0, 0, 0, 0), c(0, 0, 0, 0, 0, 0, 1, 0, 0), c(0, 0, 0, 0, 0, 1, 2, 1, 0), c(0, 0, 0, 0, 0, 0, 1, 0, 0), 0)) pm <- pivot_measure(x, y, p.1 = 1, p.2 = 2, dual.method = "epsilon-discrete") print(pm) plot(pm) # use pivot measure to calculate cost and plan pm <- transport_cost(pm) pm <- transport_df(pm) print(pm) # calculate pivot from plan pm2 <- pivot_measure(pm$df, x, y) plot(pm2)
Plot two-dimensional grids or visualize the pivot measure.
## S3 method for class 'otgrid' plot(x, num.col = 256, useRaster = TRUE, ...) ## S3 method for class 'otgridtransport' plot(x, num.col = 256, useRaster = TRUE, back.col = "lightblue", ...)
## S3 method for class 'otgrid' plot(x, num.col = 256, useRaster = TRUE, ...) ## S3 method for class 'otgridtransport' plot(x, num.col = 256, useRaster = TRUE, back.col = "lightblue", ...)
x |
an object of class |
num.col |
number of colors (grey) to use. |
useRaster |
parameter passed to the |
... |
further arguments (currently unused). |
back.col |
color of the background. |
For objects of class "otgrid"
, the grid is plotted as a greyscale image.
For objects of class "otgridtransport"
, the pivot measure is visualized as follows: the two grids of the transport are plotted
in the left upper (from
) and right lower (to
) corner. The pivot measure is in the left lower corner such that the marginals match.
No return value, called for side effects.
pivot measure pivot_measure
, grid otgrid
Calculate the optimal transport cost.
## S3 method for class 'numeric' transport_cost(x, y, wx, wy, p = 1, sorted = FALSE, threshold = 1e-15, ...) transport_cost(x, ...) ## S3 method for class 'otgridtransport' transport_cost(x, threshold = 1e-15, ...) ## S3 method for class 'otgrid' transport_cost(x, ...) ## S3 method for class 'data.frame' transport_cost(x, costm, ...)
## S3 method for class 'numeric' transport_cost(x, y, wx, wy, p = 1, sorted = FALSE, threshold = 1e-15, ...) transport_cost(x, ...) ## S3 method for class 'otgridtransport' transport_cost(x, threshold = 1e-15, ...) ## S3 method for class 'otgrid' transport_cost(x, ...) ## S3 method for class 'data.frame' transport_cost(x, costm, ...)
x |
a vector of points; a data frame with columns |
y |
second vector of points. |
wx |
weight vector of the first vector of points. |
wy |
weight vector of the second vector of points. |
p |
the power |
sorted |
logical value indicating whether or not |
threshold |
small value that indicates when a value is considered to be zero. |
... |
further arguments (for |
costm |
cost matrix of the transport |
In case of two-dimensional grids, the pivot measure is used to calculate the optimal transport cost.
For one-dimensional optimal transport, the cost function is given by . In this case, the north-west-corner
algorithm is used.
the optimal transport cost or, in case of two-dimensional case, an object of class "otgridtransport"
with element cost
that contains it.
pivot measure pivot_measure
## one-dimensional example set.seed(1) a <- 1:5 wa <- rep(1/5, 5) b <- 1:6 wb <- runif(6) wb <- wb / sum(wb) transport_cost(a, b, wa, wb, p = 1) ## two-dimensional example x <- otgrid(cbind(0:1, 1:0)) y <- otgrid(cbind(1:0, 0:1)) # first calculate pivot manually pm <- pivot_measure(x, y) pm <- transport_cost(pm) print(pm$cost) # or just pm2 <- transport_cost(x, y) print(pm2$cost) # or from transport plan and cost matrix costm <- transport_costmat(pm) tp <- transport_df(pm) print(transport_cost(tp$df, costm))
## one-dimensional example set.seed(1) a <- 1:5 wa <- rep(1/5, 5) b <- 1:6 wb <- runif(6) wb <- wb / sum(wb) transport_cost(a, b, wa, wb, p = 1) ## two-dimensional example x <- otgrid(cbind(0:1, 1:0)) y <- otgrid(cbind(1:0, 0:1)) # first calculate pivot manually pm <- pivot_measure(x, y) pm <- transport_cost(pm) print(pm$cost) # or just pm2 <- transport_cost(x, y) print(pm2$cost) # or from transport plan and cost matrix costm <- transport_costmat(pm) tp <- transport_df(pm) print(transport_cost(tp$df, costm))
Calculate the cost matrix for the optimal transport between two-dimensional grids with respect to the cost function
transport_costmat(x, ...) ## S3 method for class 'otgrid' transport_costmat(x, y, p.1 = 2, p.2 = p.1, ...) ## S3 method for class 'otgridtransport' transport_costmat(x, ...)
transport_costmat(x, ...) ## S3 method for class 'otgrid' transport_costmat(x, y, p.1 = 2, p.2 = p.1, ...) ## S3 method for class 'otgridtransport' transport_costmat(x, ...)
x |
an object of class |
... |
further arguments (currently unused). |
y |
an object of class |
p.1 |
the first power |
p.2 |
the second power |
a matrix giving the pairwise costs in column-mayor format.
x <- otgrid(cbind(0:1, 1:0)) y <- otgrid(cbind(1:0, 0:1)) transport_costmat(x, y, p.1 = 1, p.2 = 3)
x <- otgrid(cbind(0:1, 1:0)) y <- otgrid(cbind(1:0, 0:1)) transport_costmat(x, y, p.1 = 1, p.2 = 3)
Calculate the optimal transport plan.
## S3 method for class 'numeric' transport_df(x, y, threshold = 1e-15, ...) transport_df(x, ...) ## S3 method for class 'otgridtransport' transport_df(x, ...) ## S3 method for class 'otgrid' transport_df(x, ...)
## S3 method for class 'numeric' transport_df(x, y, threshold = 1e-15, ...) transport_df(x, ...) ## S3 method for class 'otgridtransport' transport_df(x, ...) ## S3 method for class 'otgrid' transport_df(x, ...)
x |
a vector of weights, an object of class |
y |
second weight vector. |
threshold |
small value that indicates when a value is considered to be zero. |
... |
further arguments (for |
In case of two-dimensional grids, the pivot measure is used to calculate the optimal transport plan.
For one-dimensional optimal transport, we assume that the optimal transport plan is the monotone plan. For example, this is the case
for costs of the form for some
.
In this case, the north-west-corner algorithm is used.
a data frame representing the optimal transport plan. It has columns from
, to
and mass
that specify
between which points of the two grids how much mass is transported.
In the two-dimensional case, a point is given by the index in column-mayor format and the data frame is actually stored in the element
df
of an object of class "otgridtransport"
.
pivot measure pivot_measure
## one-dimensional example set.seed(1) wa <- rep(1/5, 5) wb <- runif(6) wb <- wb / sum(wb) transport_df(wa, wb) ## two-dimensional example x <- otgrid(cbind(0:1, 1:0)) y <- otgrid(cbind(1:0, 0:1)) # first calculate pivot manually pm <- pivot_measure(x, y) pm <- transport_df(pm) # or just pm2 <- transport_df(x, y)
## one-dimensional example set.seed(1) wa <- rep(1/5, 5) wb <- runif(6) wb <- wb / sum(wb) transport_df(wa, wb) ## two-dimensional example x <- otgrid(cbind(0:1, 1:0)) y <- otgrid(cbind(1:0, 0:1)) # first calculate pivot manually pm <- pivot_measure(x, y) pm <- transport_df(pm) # or just pm2 <- transport_df(x, y)