Skip to contents

Given \(K\) empirical measures \(\mu_1, \mu_2, \ldots, \mu_K\) of possibly different cardinalities, wasserstein barycenter \(\mu^*\) is the solution to the following problem $$\sum_{k=1}^K \pi_k \mathcal{W}_p^p (\mu, \mu_k)$$ where \(\pi_k\)'s are relative weights of empirical measures. Here we assume either (1) support atoms in Euclidean space are given, or (2) all pairwise distances between atoms of the fixed support and empirical measures are given. The historical method used iterative Bregman projections. The current compatibility implementation is described in Details.

Usage

fbary15B(
  support,
  atoms,
  marginals = NULL,
  weights = NULL,
  lambda = 0.1,
  p = 2,
  ...
)

fbary15Bdist(
  distances,
  marginals = NULL,
  weights = NULL,
  lambda = 0.1,
  p = 2,
  ...
)

Arguments

support

an \((N\times P)\) matrix of rows being atoms for the fixed support.

atoms

a length-\(K\) list where each element is an \((N_k \times P)\) matrix of atoms.

marginals

marginal distribution for empirical measures; if NULL (default), uniform weights are set for all measures. Otherwise, it should be a length-\(K\) list where each element is a length-\(N_i\) vector of nonnegative weights that sum to 1.

weights

weights for each individual measure; if NULL (default), each measure is considered equally. Otherwise, it should be a length-\(K\) vector.

lambda

regularization parameter (default: 0.1).

p

an exponent for the order of the distance (default: 2).

...

extra parameters including

abstol

stopping criterion for iterations (default: 1e-10).

init.vec

an initial vector (default: uniform weight).

maxiter

maximum number of iterations (default: 496).

print.progress

a logical to show current iteration (default: FALSE).

distances

a length-\(K\) list where each element is an \((N\times N_k)\) pairwise distance between atoms of the fixed support and given measures.

Value

A numeric t4transport_mass vector with a diagnostics attribute and a compatible barplot method. Use return_result = TRUE to return the full fit instead.

Details

Since version 0.2.0, the historical entry points use the shared objective-checked mirror-descent solver with log-domain Sinkhorn inner solves. Their names are retained for compatibility, not as a claim that the original accelerated or Bregman-projection iteration is still used. The optimized cost is distance raised to p, plus lambda * sum(plan * (log(plan) - 1)). Locations are fixed and candidate probability masses vary. Inspect status and the reported convex objective bound; exhausting the iteration budget is not success. Additional controls are rtol, step0, max_backtrack, inner_maxiter, inner_tol, and return_result. Initial masses must be strictly positive.

References

Benamou J, Carlier G, Cuturi M, Nenna L, Peyré G (2015). “Iterative Bregman Projections for Regularized Transportation Problems.” SIAM Journal on Scientific Computing, 37(2), A1111-A1138. ISSN 1064-8275, 1095-7197, doi:10.1137/141000439 .

Examples

# \donttest{
#-------------------------------------------------------------------
#     Wasserstein Barycenter for Fixed Atoms with Two Gaussians
#
# * class 1 : samples from Gaussian with mean=(-4, -4)
# * class 2 : samples from Gaussian with mean=(+4, +4)
# * target support consists of 7 integer points from -6 to 6,
#   where ideally, weight is concentrated near 0 since it's average!
#-------------------------------------------------------------------
## GENERATE DATA
#  Empirical Measures
set.seed(100)
ndat = 20
dat1 = matrix(rnorm(ndat*2, mean=-4, sd=0.5),ncol=2)
dat2 = matrix(rnorm(ndat*2, mean=+4, sd=0.5),ncol=2) 

myatoms = list()
myatoms[[1]] = dat1
myatoms[[2]] = dat2
mydata = rbind(dat1, dat2)

#  Fixed Support
support = cbind(seq(from=-8,to=8,by=2),
                seq(from=-8,to=8,by=2))
## COMPUTE
comp1 = fbary15B(support, myatoms, lambda=0.5, maxiter=10)
#> Warning: fbary15B finished with status 'iteration_limit': maximum_iterations. Inspect attr(result, 'diagnostics').
comp2 = fbary15B(support, myatoms, lambda=1,   maxiter=10)
#> Warning: fbary15B finished with status 'iteration_limit': maximum_iterations. Inspect attr(result, 'diagnostics').
comp3 = fbary15B(support, myatoms, lambda=10,  maxiter=10)
#> Warning: fbary15B finished with status 'iteration_limit': maximum_iterations. Inspect attr(result, 'diagnostics').

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
barplot(comp1, ylim=c(0,1), main="Probability\n (lambda=0.5)")
barplot(comp2, ylim=c(0,1), main="Probability\n (lambda=1)")
barplot(comp3, ylim=c(0,1), main="Probability\n (lambda=10)")

par(opar)
# }