Package 'T4transport'

Description

digit3 contains 2000 images from the famous MNIST dataset of digit 3. Each element of the list is an image represented as an $(28\times 28)$ matrix that sums to 1. This normalization is conventional and it does not hurt its visualization via a basic 'image()' function.

Usage

data(digit3)

Format

a length-$2000$ named list "digit3" of $(28\times 28)$ matrices.

Examples

## LOAD THE DATA
data(digit3)

## SHOW A FEW
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,4), pty="s")
for (i in 1:8){
  image(digit3[[i]])
}
par(opar)

Description

digits contains 5000 images from the famous MNIST dataset of all digits, consisting of 500 images per digit class from 0 to 9. Each digit image is represented as an $(28\times 28)$ matrix that sums to 1. This normalization is conventional and it does not hurt its visualization via a basic 'image()' function.

Usage

data(digits)

Format

a named list "digits" containing

image

length-5000 list of $(28\times 28)$ image matrices.

label

length-5000 vector of class labels from 0 to 9.

Examples

## LOAD THE DATA
data(digits)

## SHOW A FEW
#  Select 9 random images
subimgs = digits$image[sample(1:5000, 9)]

opar <- par(no.readonly=TRUE)
par(mfrow=c(3,3), pty="s")
for (i in 1:9){
  image(subimgs[[i]])
}
par(opar)

Description

Given a collection of empirical cumulative distribution functions $F^i (x)$ for $i=1,\ldots,N$, compute the Wasserstein barycenter of order 2. This is obtained by taking a weighted average on a set of corresponding quantile functions.

Usage

ecdfbary(ecdfs, weights = NULL, ...)

Arguments

Value

an "ecdf" object of the Wasserstein barycenter.

Examples

#----------------------------------------------------------------------
#                         Two Gaussians
#
# Two Gaussian distributions are parametrized as follows.
# Type 1 : (mean, var) = (-4, 1/4)
# Type 2 : (mean, var) = (+4, 1/4)
#----------------------------------------------------------------------
# GENERATE ECDFs
ecdf_list = list()
ecdf_list[[1]] = stats::ecdf(stats::rnorm(200, mean=-4, sd=0.5))
ecdf_list[[2]] = stats::ecdf(stats::rnorm(200, mean=+4, sd=0.5))

# COMPUTE THE BARYCENTER OF EQUAL WEIGHTS
emean = ecdfbary(ecdf_list)

# QUANTITIES FOR PLOTTING
x_grid  = seq(from=-8, to=8, length.out=100)
y_type1 = ecdf_list[[1]](x_grid)
y_type2 = ecdf_list[[2]](x_grid)
y_bary  = emean(x_grid)

# VISUALIZE
opar <- par(no.readonly=TRUE)
plot(x_grid, y_bary, lwd=3, col="red", type="l",
     main="Barycenter", xlab="x", ylab="Fn(x)")
lines(x_grid, y_type1, col="gray50", lty=3)
lines(x_grid, y_type2, col="gray50", lty=3)
par(opar)

Description

Given a collection of empirical cumulative distribution functions $F^i (x)$ for $i=1,\ldots,N$, compute the Wasserstein median. This is obtained by a functional variant of the Weiszfeld algorithm on a set of quantile functions.

Usage

ecdfmed(ecdfs, weights = NULL, ...)

Arguments

Value

an "ecdf" object of the Wasserstein median.

Examples

#----------------------------------------------------------------------
#                         Tree Gaussians
#
# Three Gaussian distributions are parametrized as follows.
# Type 1 : (mean, sd) = (-4, 1)
# Type 2 : (mean, sd) = ( 0, 1/5)
# Type 3 : (mean, sd) = (+6, 1/2)
#----------------------------------------------------------------------
# GENERATE ECDFs
ecdf_list = list()
ecdf_list[[1]] = stats::ecdf(stats::rnorm(200, mean=-4, sd=1))
ecdf_list[[2]] = stats::ecdf(stats::rnorm(200, mean=+4, sd=0.2))
ecdf_list[[3]] = stats::ecdf(stats::rnorm(200, mean=+6, sd=0.5))

# COMPUTE THE MEDIAN
emeds = ecdfmed(ecdf_list)

# COMPUTE THE BARYCENTER
emean = ecdfbary(ecdf_list)

# QUANTITIES FOR PLOTTING
x_grid  = seq(from=-8, to=10, length.out=500)
y_type1 = ecdf_list[[1]](x_grid)
y_type2 = ecdf_list[[2]](x_grid)
y_type3 = ecdf_list[[3]](x_grid)

y_bary = emean(x_grid)
y_meds = emeds(x_grid)

# VISUALIZE
opar <- par(no.readonly=TRUE)
plot(x_grid, y_bary, lwd=3, col="orange", type="l",
     main="Wasserstein Median & Barycenter", 
     xlab="x", ylab="Fn(x)", lty=2)
lines(x_grid, y_meds, lwd=3, col="blue", lty=2)
lines(x_grid, y_type1, col="gray50", lty=3)
lines(x_grid, y_type2, col="gray50", lty=3)
lines(x_grid, y_type3, col="gray50", lty=3)
legend("topleft", legend=c("Median","Barycenter"),
        lwd=3, lty=2, col=c("blue","orange"))
par(opar)

Description

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. Algorithmically, it is a subgradient method where the each subgradient is approximated using the entropic regularization.

Usage

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

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

Arguments

Value

a length-$N$ vector of probability vector.

References

Cuturi M, Doucet A (2014-06-22/2014-06-24). “Fast Computation of Wasserstein Barycenters.” In Xing EP, Jebara T (eds.), Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, 685–693.

Examples

#-------------------------------------------------------------------
#     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 = 100
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 = fbary14C(support, myatoms, lambda=0.5, maxiter=10)
comp2 = fbary14C(support, myatoms, lambda=1,   maxiter=10)
comp3 = fbary14C(support, myatoms, lambda=10,  maxiter=10)

## 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)

Description

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. Authors proposed iterative Bregman projections in conjunction with entropic regularization.

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

Value

a length-$N$ vector of probability vector.

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

#-------------------------------------------------------------------
#     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 = 500
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)
comp2 = fbary15B(support, myatoms, lambda=1,   maxiter=10)
comp3 = fbary15B(support, myatoms, lambda=10,  maxiter=10)

## 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)

Description

Given a point cloud $X \in \mathbf{R}^{N \times P}$, this function constructs a fully connected weighted graph using an RBF (Gaussian) kernel with bandwidth chosen by the median heuristic, forms the unnormalized graph Laplacian, and returns the corresponding Fiedler vector, which is the eigenvector associated to the second smallest eigenvalue of the Laplacian.

Usage

fiedler(X, normalize = TRUE)

Arguments

Value

A numeric vector of length $N$ containing the Fiedler values associated with each point in the input point cloud. If normalize = TRUE, the entries are in the interval $[0,1]$.

Examples

#-------------------------------------------------------------------
#                           Description
#
# Use 'iris' dataset to compute fiedler vector. 
# The dataset is visualized in R^2 using PCA
#-------------------------------------------------------------------
# load dataset
X = as.matrix(iris[,1:4])

# PCA preprocessing
X2d = X%*%eigen(cov(X))$vectors[,1:2]

# compute fiedler vector
fied_vec = fiedler(X2d, normalize=TRUE)

# plot 
opar <- par(no.readonly=TRUE)
plot(X2d, col=rainbow(150)[as.numeric(cut(fied_vec, breaks=150))], 
     pch=19, xlab="PC 1", ylab="PC 2",
     main="Fiedler vector on Iris dataset (PCA-reduced)")
par(opar)

Description

Given a collection of Gaussian distributions $\mathcal{N}(\mu_i, \sigma_i^2)$ for $i=1,\ldots,n$, compute the Wasserstein barycenter of order 2. For the barycenter computation of variance components, we use a fixed-point algorithm by Álvarez-Esteban et al. (2016).

Usage

gaussbary1d(means, vars, weights = NULL, ...)

Arguments

Value

a named list containing

mean

mean of the estimated barycenter distribution.

var

variance of the estimated barycenter distribution.

References

Álvarez-Esteban PC, del Barrio E, Cuesta-Albertos JA, Matrán C (2016). “A Fixed-Point Approach to Barycenters in Wasserstein Space.” Journal of Mathematical Analysis and Applications, 441(2), 744–762. ISSN 0022247X. doi:10.1016/j.jmaa.2016.04.045.

See Also

[T4transport::gaussbarypd()] for multivariate case.

Examples

#----------------------------------------------------------------------
#                         Two Gaussians
#
# Two Gaussian distributions are parametrized as follows.
# Type 1 : (mean, var) = (-4, 1/4)
# Type 2 : (mean, var) = (+4, 1/4)
#----------------------------------------------------------------------
# GENERATE PARAMETERS
par_mean = c(-4, 4)
par_vars = c(0.25, 0.25)

# COMPUTE THE BARYCENTER OF EQUAL WEIGHTS
gmean = gaussbary1d(par_mean, par_vars)

# QUANTITIES FOR PLOTTING
x_grid  = seq(from=-6, to=6, length.out=200)
y_dist1 = stats::dnorm(x_grid, mean=-4, sd=0.5)
y_dist2 = stats::dnorm(x_grid, mean=+4, sd=0.5)
y_gmean = stats::dnorm(x_grid, mean=gmean$mean, sd=sqrt(gmean$var)) 

# VISUALIZE
opar <- par(no.readonly=TRUE)
plot(x_grid, y_gmean, lwd=2, col="red", type="l",
     main="Barycenter", xlab="x", ylab="density")
lines(x_grid, y_dist1)
lines(x_grid, y_dist2)
par(opar)

Description

Given a collection of $n$-dimensional Gaussian distributions $N(\mu_i, \Sigma_i)$ for $i=1,\ldots,n$, compute the Wasserstein barycenter of order 2. For the barycenter computation of variance components, we use a fixed-point algorithm by Álvarez-Esteban et al. (2016).

Usage

gaussbarypd(means, vars, weights = NULL, ...)

Arguments

Value

a named list containing

mean

a length-$p$ vector for mean of the estimated barycenter distribution.

var

a $(p\times p)$ matrix for variance of the estimated barycenter distribution.

References

Álvarez-Esteban PC, del Barrio E, Cuesta-Albertos JA, Matrán C (2016). “A Fixed-Point Approach to Barycenters in Wasserstein Space.” Journal of Mathematical Analysis and Applications, 441(2), 744–762. ISSN 0022247X. doi:10.1016/j.jmaa.2016.04.045.

See Also

[T4transport::gaussbary1d()] for univariate case.

Examples

#----------------------------------------------------------------------
#                         Two Gaussians in R^2
#----------------------------------------------------------------------
# GENERATE PARAMETERS
# means
par_mean = rbind(c(-4,0), c(4,0))

# covariances
par_vars = array(0,c(2,2,2))
par_vars[,,1] = cbind(c(4,-2),c(-2,4))
par_vars[,,2] = cbind(c(4,+2),c(+2,4))

# COMPUTE THE BARYCENTER OF EQUAL WEIGHTS
gmean = gaussbarypd(par_mean, par_vars)

# GET COORDINATES FOR DRAWING
pt_type1 = gaussvis2d(par_mean[1,], par_vars[,,1])
pt_type2 = gaussvis2d(par_mean[2,], par_vars[,,2])
pt_gmean = gaussvis2d(gmean$mean, gmean$var)

# VISUALIZE
opar <- par(no.readonly=TRUE)
plot(pt_gmean, lwd=2, col="red", type="l",
     main="Barycenter", xlab="", ylab="", 
     xlim=c(-6,6))
lines(pt_type1)
lines(pt_type2)
par(opar)

Description

Given a collection of Gaussian distributions $\mathcal{N}(\mu_i, \sigma_i^2)$ for $i=1,\ldots,n$, compute the Wasserstein median.

Usage

gaussmed1d(means, vars, weights = NULL, ...)

Arguments

Value

a named list containing

mean

mean of the estimated median distribution.

var

variance of the estimated median distribution.

References

You K, Shung D, Giuffrè M (2025). “On the Wasserstein Median of Probability Measures.” Journal of Computational and Graphical Statistics, 34(1), 253-266. ISSN 1061-8600, 1537-2715.

See Also

[T4transport::gaussmedpd()] for multivariate case.

Examples

#----------------------------------------------------------------------
#                         Tree Gaussians
#
# Three Gaussian distributions are parametrized as follows.
# Type 1 : (mean, sd) = (-4, 1)
# Type 2 : (mean, sd) = ( 0, 1/5)
# Type 3 : (mean, sd) = (+6, 1/2)
#----------------------------------------------------------------------
# GENERATE PARAMETERS
par_mean = c(-4, 0, +6)
par_vars = c(1, 0.04, 0.25)

# COMPUTE THE WASSERSTEIN MEDIAN
gmeds = gaussmed1d(par_mean, par_vars)

# COMPUTE THE BARYCENTER 
gmean = gaussbary1d(par_mean, par_vars)

# QUANTITIES FOR PLOTTING
x_grid  = seq(from=-6, to=8, length.out=1000)
y_dist1 = stats::dnorm(x_grid, mean=par_mean[1], sd=sqrt(par_vars[1]))
y_dist2 = stats::dnorm(x_grid, mean=par_mean[2], sd=sqrt(par_vars[2]))
y_dist3 = stats::dnorm(x_grid, mean=par_mean[3], sd=sqrt(par_vars[3]))

y_gmean = stats::dnorm(x_grid, mean=gmean$mean, sd=sqrt(gmean$var)) 
y_gmeds = stats::dnorm(x_grid, mean=gmeds$mean, sd=sqrt(gmeds$var))

# VISUALIZE
opar <- par(no.readonly=TRUE)
plot(x_grid, y_gmeds, lwd=3, col="red", type="l",
     main="Three Gaussians", xlab="x", ylab="density", 
     xlim=range(x_grid), ylim=c(0,2.5))
lines(x_grid, y_gmean, lwd=3, col="blue")
lines(x_grid, y_dist1, lwd=1.5, lty=2)
lines(x_grid, y_dist2, lwd=1.5, lty=2)
lines(x_grid, y_dist3, lwd=1.5, lty=2)
legend("topleft", legend=c("Median","Barycenter"),
       col=c("red","blue"), lwd=c(3,3), lty=c(1,2))
par(opar)

Description

Given a collection of $p$-dimensional Gaussian distributions $N(\mu_i, \Sigma_i)$ for $i=1,\ldots,n$, compute the Wasserstein median.

Usage

gaussmedpd(means, vars, weights = NULL, ...)

Arguments

Value

a named list containing

mean

a length-$p$ vector for mean of the estimated median distribution.

var

a $(p\times p)$ matrix for variance of the estimated median distribution.

References

You K, Shung D, Giuffrè M (2025). “On the Wasserstein Median of Probability Measures.” Journal of Computational and Graphical Statistics, 34(1), 253-266. ISSN 1061-8600, 1537-2715.

See Also

[T4transport::gaussmed1d()] for univariate case.

Examples

#----------------------------------------------------------------------
#                         Three Gaussians in R^2
#----------------------------------------------------------------------
# GENERATE PARAMETERS
# means
par_mean = rbind(c(-4,0), c(0,0), c(5,-1))

# covariances
par_vars = array(0,c(2,2,3))
par_vars[,,1] = cbind(c(2,-1),c(-1,2))
par_vars[,,2] = cbind(c(4,+1),c(+1,4))
par_vars[,,3] = diag(c(4,1))

# COMPUTE THE MEDIAN
gmeds = gaussmedpd(par_mean, par_vars)

# COMPUTE THE BARYCENTER 
gmean = gaussbarypd(par_mean, par_vars)

# GET COORDINATES FOR DRAWING
pt_type1 = gaussvis2d(par_mean[1,], par_vars[,,1])
pt_type2 = gaussvis2d(par_mean[2,], par_vars[,,2])
pt_type3 = gaussvis2d(par_mean[3,], par_vars[,,3])
pt_gmean = gaussvis2d(gmean$mean, gmean$var)
pt_gmeds = gaussvis2d(gmeds$mean, gmeds$var)

# VISUALIZE
opar <- par(no.readonly=TRUE)
plot(pt_gmean, lwd=2, col="red", type="l",
     main="Three Gaussians", xlab="", ylab="", 
     xlim=c(-6,8), ylim=c(-2.5,2.5))
lines(pt_gmeds, lwd=2, col="blue")
lines(pt_type1, lty=2, lwd=5)
lines(pt_type2, lty=2, lwd=5)
lines(pt_type3, lty=2, lwd=5)
abline(h=0, col="grey80", lty=3)
abline(v=0, col="grey80", lty=3)
legend("topright", legend=c("Median","Barycenter"),
       lwd=2, lty=1, col=c("blue","red"))
par(opar)

Description

This function samples points along the contour of an ellipse represented by mean and variance parameters for a 2-dimensional Gaussian distribution to help ease manipulating visualization of the specified distribution. For example, you can directly use a basic plot() function directly for drawing.

Usage

gaussvis2d(mean, var, n = 500)

Arguments

Value

an $(n\times 2)$ matrix.

Examples

#----------------------------------------------------------------------
#                        Three Gaussians in R^2
#----------------------------------------------------------------------
# MEAN PARAMETERS
loc1 = c(-3,0)
loc2 = c(0,5)
loc3 = c(3,0)

# COVARIANCE PARAMETERS
var1 = cbind(c(4,-2),c(-2,4))
var2 = diag(c(9,1))
var3 = cbind(c(4,2),c(2,4))

# GENERATE POINTS
visA = gaussvis2d(loc1, var1)
visB = gaussvis2d(loc2, var2)
visC = gaussvis2d(loc3, var3)

# VISUALIZE
opar <- par(no.readonly=TRUE)
plot(visA[,1], visA[,2], type="l", xlim=c(-5,5), ylim=c(-2,9),
     lwd=3, col="red", main="3 Gaussian Distributions")
lines(visB[,1], visB[,2], lwd=3, col="blue")
lines(visC[,1], visC[,2], lwd=3, col="orange")
legend("top", legend=c("Type 1","Type 2","Type 3"),
       lwd=3, col=c("red","blue","orange"), horiz=TRUE)
par(opar)

Description

Computes the Gromov–Wasserstein (GW) barycenter of a collection of metric measure spaces. Given a list of distance matrices ${D^{(k)}}{k=1}^K$ and their corresponding marginal distributions, the function estimates a synthetic metric space whose intrinsic geometry best represents the input collection under the GW criterion.

The GW barycenter is defined as the minimizer of a multi-measure Gromov–Wasserstein objective, where each dataset contributes according to a user-specified barycentric weight. Since the problem is jointly non-convex in the barycenter metric and the coupling matrices, the algorithm proceeds through an outer–inner iterative procedure.

Usage

gwbary(distances, marginals = NULL, weights = NULL, num_support = 100, ...)

Arguments

Value

A named list containing

dist

an object of class dist representing the GW barycenter.

weight

a length-$M$ vector of barycenter weights with all entries being $1/M$.

Examples

## Not run: 
#-------------------------------------------------------------------
#                           Description
#
# GW barycenter computation is quite expensive. In this example, 
# we draw a small set of empirical measures from the digit '3'
# images and compute their GW barycenter with a small number of 
# support points. The attained barycenter distance matrix is then 
# passed onto the classical MDS algorithm for visualization.
#-------------------------------------------------------------------
## GENERATE DATA
data(digits)
data_D = vector("list", length=5)
data_W = vector("list", length=5)
for (i in 1:5){
  img_now = img2measure(digits3[[i]])
  data_D[[i]] = stats::dist(img_now$support)
  data_W[[i]] = as.vector(img_now$weight)
}

## COMPUTE
bary_dist <- gwbary(data_D, marginals=data_W, num_support=100)
bary_cmd2 <- stats::cmdscale(bary_dist$dist, k=2)

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(pty="s")
plot(bary_cmd2, main="GW Barycenter Embedding",
     xaxt="n", yaxt="n", pch=19, xlab="", ylab="")
par(opar)

## End(Not run)

Description

Computes the Gromov-Wasserstein (GW) distance between two metric measure spaces. Given two distance matrices $D_X$ and $D_Y$ along with their respective marginal distributions, the function solves the GW optimization problem to obtain both the distance value and an associated optimal transport plan.

The GW distance provides a way to compare datasets that may not lie in the same ambient space by focusing on the intrinsic geometric structure encoded in the pairwise distances. This implementation supports multiple optimization schemes, including majorization–minimization (MM), proximal gradient (PG), and Frank–Wolfe (FW).

Usage

gwdist(Dx, Dy, wx = NULL, wy = NULL, ...)

Arguments

Value

a named list containing

distance

the computed GW distance value.

plan

an $(M\times N)$ nonnegative matrix for the optimal transport plan.

References

Mémoli F (2011). “Gromov–Wasserstein Distances and the Metric Approach to Object Matching.” Foundations of Computational Mathematics, 11(4), 417–487. ISSN 1615-3375, 1615-3383.

Examples

## Not run: 
#-------------------------------------------------------------------
#                           Description
#
# * class 1 : iris dataset (columns 1-4) with perturbations
# * class 2 : class 1 rotated randomly in R^4
# * class 3 : samples from N((0,0), I)
#
#  We draw 10 empirical measures from each and compare 
#  the regular Wasserstein and GW distance. It is expected that 
#  the GW distance between class 1 and class 2 is negligible, 
#  while the regular Wasserstein distance is large. For simplicity, 
#  limit the cardinalities to 20.
#-------------------------------------------------------------------
## GENERATE DATA
set.seed(10)

#  prepare empty lists
inputs = vector("list", length=30)

#  generate class 1 and 2
iris_mat = as.matrix(iris[sample(1:150,20),1:4])
for (i in 1:10){
  inputs[[i]] = iris_mat + matrix(rnorm(20*4), ncol=4)
  inputs[[i+10]] = inputs[[i]]%*%qr.Q(qr(matrix(runif(16), ncol=4)))
}
#  generate class 3
for (j in 21:30){
  inputs[[j]] = matrix(rnorm(20*4), ncol=4)
}

## COMPUTE
#  empty arrays
dist_RW = array(0, c(30, 30))
dist_GW = array(0, c(30, 30))

#  compute pairwise distances
for (i in 1:29){
  X <- inputs[[i]]
  Dx <- stats::dist(X)
  for (j in (i+1):30){
  Y <- inputs[[j]]
  Dy <- stats::dist(Y)
  dist_RW[i,j] <- dist_RW[j,i] <- wasserstein(X, Y)$distance
  dist_GW[i,j] <- dist_GW[j,i] <- gwdist(Dx, Dy)$distance
  }
}

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2), pty="s")
image(dist_RW, xaxt="n", yaxt="n", main="Regular Wasserstein distance")
image(dist_GW, xaxt="n", yaxt="n", main="Gromov-Wasserstein distance")
par(opar)

## End(Not run)

Description

Given multiple histograms represented as "histogram" S3 objects, compute their 2-Wasserstein barycenter using the exact 1D quantile characterization. All input histograms must have identical breaks.

Usage

histbary(hists, weights = NULL, L = 2000L)

Arguments

Value

a "histogram" object representing the Wasserstein barycenter.

Examples

#----------------------------------------------------------------------
#                      Binned from Two Gaussians
#
# EXAMPLE : Very Small Example for CRAN; just showing how to use it!
#----------------------------------------------------------------------
# GENERATE FROM TWO GAUSSIANS WITH DIFFERENT MEANS
set.seed(100)
x  = stats::rnorm(1000, mean=-4, sd=0.5)
y  = stats::rnorm(1000, mean=+4, sd=0.5)
bk = seq(from=-10, to=10, length.out=20)

# HISTOGRAMS WITH COMMON BREAKS
histxy = list()
histxy[[1]] = hist(x, breaks=bk, plot=FALSE)
histxy[[2]] = hist(y, breaks=bk, plot=FALSE)

# COMPUTE
hh = histbary(histxy)

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2), pty="s")
barplot(histxy[[1]]$density, col=rgb(0,0,1,1/4), 
        ylim=c(0, 0.75), main="Two Histograms")
barplot(histxy[[2]]$density, col=rgb(1,0,0,1/4), 
        ylim=c(0, 0.75), add=TRUE)
barplot(hh$density, main="Barycenter",
        ylim=c(0, 0.75))
par(opar)

Description

Given multiple histograms represented as "histogram" S3 objects, compute Wasserstein barycenter. We need one requirement that all histograms in an input list hists must have same breaks. See the example on how to construct a histogram on predefined breaks/bins.

Usage

histbary14C(hists, p = 2, weights = NULL, lambda = NULL, ...)

Arguments

Value

a "histogram" object representing the Wasserstein barycenter.

References

Cuturi M, Doucet A (2014-06-22/2014-06-24). “Fast Computation of Wasserstein Barycenters.” In Xing EP, Jebara T (eds.), Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, 685–693.

See Also

fbary14C

Examples

#----------------------------------------------------------------------
#                      Binned from Two Gaussians
#
# EXAMPLE : Very Small Example for CRAN; just showing how to use it!
#----------------------------------------------------------------------
# GENERATE FROM TWO GAUSSIANS WITH DIFFERENT MEANS
set.seed(100)
x  = stats::rnorm(1000, mean=-4, sd=0.5)
y  = stats::rnorm(1000, mean=+4, sd=0.5)
bk = seq(from=-10, to=10, length.out=20)

# HISTOGRAMS WITH COMMON BREAKS
histxy = list()
histxy[[1]] = hist(x, breaks=bk, plot=FALSE)
histxy[[2]] = hist(y, breaks=bk, plot=FALSE)

# COMPUTE
hh = histbary14C(histxy, maxiter=5)

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2))
barplot(histxy[[1]]$density, col=rgb(0,0,1,1/4), 
        ylim=c(0, 0.75), main="Two Histograms")
barplot(histxy[[2]]$density, col=rgb(1,0,0,1/4), 
        ylim=c(0, 0.75), add=TRUE)
barplot(hh$density, main="Barycenter",
        ylim=c(0, 0.75))
par(opar)

Description

Given multiple histograms represented as "histogram" S3 objects, compute Wasserstein barycenter. We need one requirement that all histograms in an input list hists must have same breaks. See the example on how to construct a histogram on predefined breaks/bins.

Usage

histbary15B(hists, p = 2, weights = NULL, lambda = NULL, ...)

Arguments

Value

a "histogram" object of barycenter.

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.

See Also

fbary15B

Examples

#----------------------------------------------------------------------
#                      Binned from Two Gaussians
#
# EXAMPLE : Very Small Example for CRAN; just showing how to use it!
#----------------------------------------------------------------------
# GENERATE FROM TWO GAUSSIANS WITH DIFFERENT MEANS
set.seed(100)
x  = stats::rnorm(1000, mean=-4, sd=0.5)
y  = stats::rnorm(1000, mean=+4, sd=0.5)
bk = seq(from=-10, to=10, length.out=20)

# HISTOGRAMS WITH COMMON BREAKS
histxy = list()
histxy[[1]] = hist(x, breaks=bk, plot=FALSE)
histxy[[2]] = hist(y, breaks=bk, plot=FALSE)

# COMPUTE
hh = histbary15B(histxy, maxiter=5)

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2))
barplot(histxy[[1]]$density, col=rgb(0,0,1,1/4), 
        ylim=c(0, 0.75), main="Two Histograms")
barplot(histxy[[2]]$density, col=rgb(1,0,0,1/4), 
        ylim=c(0, 0.75), add=TRUE)
barplot(hh$density, main="Barycenter",
        ylim=c(0, 0.75))
par(opar)

Description

Compute the $p$-Wasserstein distance between two 1D histograms that share the same binning, i.e., same breaks. The histograms are treated as discrete probability measures supported at bin midpoints with masses given by normalized counts. Uses the exact 1D monotone OT algorithm, not LP nor entropic regularization.

Usage

histdist(hist1, hist2, p = 2)

Arguments

Value

a named list containing

distance

$\mathcal{W}_p$ distance value.

Examples

#----------------------------------------------------------------------
#                      Binned from Gaussian and Uniform
#
# Create two types of histograms with the same binning. One is from 
# the standard normal and the other from uniform distribution in [-5,5].
#----------------------------------------------------------------------
# GENERATE 20 HISTOGRAMS
set.seed(100)
hist20 = list()
bk = seq(from=-10, to=10, length.out=20) # common breaks
for (i in 1:10){
  hist20[[i]] = hist(stats::rnorm(100), breaks=bk, plot=FALSE)
  hist20[[i+10]] = hist(stats::runif(100, min=-5, max=5), breaks=bk, plot=FALSE)
}

# COMPUTE THE PAIRWISE DISTANCE
pdmat = array(0,c(20,20))
for (i in 1:19){
  for (j in (i+1):20){
    pdmat[i,j] = histdist(hist20[[i]], hist20[[j]], p=2)$distance
    pdmat[j,i] = pdmat[i,j]
  }
}

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(pty="s")
image(pdmat, axes=FALSE, main="Pairwise 2-Wasserstein Distance between Histograms")
par(opar)

Description

Given two histograms represented as "histogram" S3 objects with identical breaks, compute interpolated histograms along the 2-Wasserstein geodesic connecting them. In 1D, this is achieved by linear interpolation of quantile functions (displacement interpolation).

Usage

histinterp(hist1, hist2, t = 0.5, L = 2000L)

Arguments

Value

If length(t) == 1, a single "histogram" object representing the interpolated distribution at time t. If length(t) > 1, a length-length(t) list of "histogram" objects.

Examples

#----------------------------------------------------------------------
#                      Interpolating Two Gaussians
#
# The source histogram is created from N(-5,1/4).
# The target histogram is created from N(+5,4)
#----------------------------------------------------------------------
# SETTING
set.seed(123)
x_source = rnorm(1000, mean=-5, sd=1/2)
x_target = rnorm(1000, mean=+5, sd=2)

# BUILD HISTOGRAMS WITH COMMON BREAKS
bk = seq(from=-8, to=12, by=2)
h1 = hist(x_source, breaks=bk, plot=FALSE)
h2 = hist(x_target, breaks=bk, plot=FALSE)

# INTERPOLATE WITH 5 GRID POINTS
h_path <- histinterp(h1, h2, t = seq(0, 1, length.out = 8))

# VISUALIZE
y_slim <- c(0, max(h1$density, h2$density)) # shared y-limit
xt     <- round(h1$mids, 1) # x-ticks

opar <- par(no.readonly = TRUE)
par(mfrow = c(2,4), pty = "s")
for (i in 1:8){
  if (i < 2){
    barplot(h_path[[i]]$density,  names.arg=xt, ylim=y_slim,
            main="Source", col=rgb(0,0,1,1/4))
  } else if (i > 7){
    barplot(h_path[[i]]$density,  names.arg=xt, ylim=y_slim,
            main="Target", col=rgb(1,0,0,1/4))
  } else {
    barplot(h_path[[i]]$density,  names.arg=xt, ylim=y_slim, 
            col="gray90", main=sprintf("t = %.3f", (i-1)/7))
  }
}
par(opar)

Description

Given multiple histograms represented as "histogram" S3 objects with common breaks, compute their Fréchet (geometric) median under the 2-Wasserstein distance. In 1D, this is implemented by mapping histograms to their quantile functions and running a Weiszfeld-type algorithm for the geometric median in the Hilbert space of quantile functions.

Usage

histmed(hists, weights = NULL, L = 2000L, ...)

Arguments

Value

a "histogram" object representing the Wasserstein median.

Examples

#----------------------------------------------------------------------
#                      Binned from Two Gaussians
#
# Generate 12 histograms from N(-4,1/4) and 8 from N(4,1/4)
#----------------------------------------------------------------------
# COMMON SETTING
set.seed(100)
bk = seq(from=-10, to=10, length.out=20)
n_signal = 12

# GENERATE HISTOGRAMS WITH COMMON BREAKS
hist_all = list()
for (i in 1:n_signal){
  hist_all[[i]] = hist(stats::rnorm(200, mean=-4, sd=0.5), breaks=bk)
}
for (j in (n_signal+1):20){
  hist_all[[j]] = hist(stats::rnorm(200, mean=+4, sd=0.5), breaks=bk)
}

# COMPUTE THE BARYCENTER AND THE MEDIAN 
h_bary = histbary(hist_all)
h_med  = histmed(hist_all)

# VISUALIZE
xt   <- round(h_med$mids, 1) 
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
barplot(hist_all[[1]]$density, col=rgb(0,0,1,1/4), 
        ylim=c(0, 0.75), main="Two Types", names.arg=xt)
barplot(hist_all[[20]]$density, col=rgb(1,0,0,1/4), 
        ylim=c(0, 0.75), add=TRUE)
barplot(h_med$density,  names.arg=xt, main="Median", ylim=c(0, 0.75))
barplot(h_bary$density, names.arg=xt, main="Barycenter", ylim=c(0, 0.75))
par(opar)

Description

Using exact balanced optimal transport as a subroutine, imagebary computes an unregularized 2-Wasserstein barycenter image $X^\star$ from multiple input images $X_1,\ldots,X_N$. Unlike the other image barycenter routines, this function does not use entropic regularization. Instead, it solves the barycenter problem with a robust first-order method based on mirror descent on the probability simplex.

Usage

imagebary(images, p = 2, weights = NULL, C = NULL, ...)

Arguments

Details

The algorithm treats each image as a discrete probability distribution on a common $(m\times n)$ grid. At each iteration, it computes exact OT dual potentials $u_i$ between the current barycenter iterate and each input image via util_dual_emd_C. These dual potentials form a valid subgradient of the barycenter objective, and a KL-mirror descent step produces a strictly positive update of the barycenter weights. For numerical stability, the implementation includes (i) centering of dual potentials (shift invariance), (ii) gradient clipping, (iii) log-domain normalization, and (iv) optional smoothing/backtracking safeguards to avoid infeasible OT calls.

Value

an $(m\times n)$ matrix of the barycentric image.

Examples

## Not run: 
#----------------------------------------------------------------------
#                       MNIST Data with Digit 3
#
# small example to compare the un- and regularized problem solutions
# choose only 10 images and run for 20 iterations with default penalties
#----------------------------------------------------------------------
# LOAD DATA
set.seed(11)
data(digit3)
dat_small = digit3[sample(1:2000, 10)]

# RUN
run_exact = imagebary(dat_small, maxiter=20)
run_reg14 = imagebary14C(dat_small, maxiter=20)
run_reg15 = imagebary15B(dat_small, maxiter=20)

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
image(run_exact, axes=FALSE, main="Unregularized")
image(run_reg14, axes=FALSE, main="Cuturi & Doucet (2014)")
image(run_reg15, axes=FALSE, main="Benamou et al. (2015)")
par(opar)

## End(Not run)

Description

Using entropic regularization for Wasserstein barycenter computation, imagebary14C finds a barycentric image $X^*$ given multiple images $X_1,X_2,\ldots,X_N$. Please note the followings; (1) we only take a matrix as an image so please make it grayscale if not, (2) all images should be of same size - no resizing is performed.

Usage

imagebary14C(images, p = 2, weights = NULL, lambda = NULL, ...)

Arguments

Value

an $(m\times n)$ matrix of the barycentric image.

References

Cuturi M, Doucet A (2014-06-22/2014-06-24). “Fast Computation of Wasserstein Barycenters.” In Xing EP, Jebara T (eds.), Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, 685–693.

See Also

fbary14C

Examples

## Not run: 
#----------------------------------------------------------------------
#                       MNIST Data with Digit 3
#
# EXAMPLE 1 : Very Small  Example for CRAN; just showing how to use it!
# EXAMPLE 2 : Medium-size Example for Evolution of Output
#----------------------------------------------------------------------
# EXAMPLE 1
data(digit3)
datsmall = digit3[1:2]
outsmall = imagebary14C(datsmall, maxiter=3)

# EXAMPLE 2 : Barycenter of 100 Images
# RANDOMLY SELECT THE IMAGES
data(digit3)
dat2 = digit3[sample(1:2000, 100)]  # select 100 images

# RUN SEQUENTIALLY
run10 = imagebary14C(dat2, maxiter=10)                   # first 10 iterations
run20 = imagebary14C(dat2, maxiter=10, init.image=run10) # run 40 more
run50 = imagebary14C(dat2, maxiter=30, init.image=run20) # run 50 more

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,3), pty="s")
image(dat2[[sample(100,1)]], axes=FALSE, main="a random image")
image(dat2[[sample(100,1)]], axes=FALSE, main="a random image")
image(dat2[[sample(100,1)]], axes=FALSE, main="a random image")
image(run10, axes=FALSE, main="barycenter after 10 iter")
image(run20, axes=FALSE, main="barycenter after 20 iter")
image(run50, axes=FALSE, main="barycenter after 50 iter")
par(opar)

## End(Not run)

Description

Using entropic regularization for Wasserstein barycenter computation, imagebary15B finds a barycentric image $X^*$ given multiple images $X_1,X_2,\ldots,X_N$. Please note the followings; (1) we only take a matrix as an image so please make it grayscale if not, (2) all images should be of same size - no resizing is performed.

Usage

imagebary15B(images, p = 2, weights = NULL, lambda = NULL, ...)

Arguments

Value

an $(m\times n)$ matrix of the barycentric image.

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.

See Also

fbary15B

Examples

#----------------------------------------------------------------------
#                       MNIST Data with Digit 3
#
# EXAMPLE 1 : Very Small  Example for CRAN; just showing how to use it!
# EXAMPLE 2 : Medium-size Example for Evolution of Output
#----------------------------------------------------------------------
# EXAMPLE 1
data(digit3)
datsmall = digit3[1:2]
outsmall = imagebary15B(datsmall, maxiter=3)

## Not run: 
# EXAMPLE 2 : Barycenter of 100 Images
# RANDOMLY SELECT THE IMAGES
data(digit3)
dat2 = digit3[sample(1:2000, 100)]  # select 100 images

# RUN SEQUENTIALLY
run05 = imagebary15B(dat2, maxiter=5)                    # first 5 iterations
run10 = imagebary15B(dat2, maxiter=5,  init.image=run05) # run 5 more
run50 = imagebary15B(dat2, maxiter=40, init.image=run10) # run 40 more

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,3), pty="s")
image(dat2[[sample(100,1)]], axes=FALSE, main="a random image")
image(dat2[[sample(100,1)]], axes=FALSE, main="a random image")
image(dat2[[sample(100,1)]], axes=FALSE, main="a random image")
image(run05, axes=FALSE, main="barycenter after 05 iter")
image(run10, axes=FALSE, main="barycenter after 10 iter")
image(run50, axes=FALSE, main="barycenter after 50 iter")
par(opar)

## End(Not run)

Description

Given two grayscale images represented as numeric matrices, compute their Wasserstein distance using an exact balanced optimal transport solver. Each image is interpreted as a discrete probability distribution on a common $(m\times n)$ grid. The ground cost is defined using the Euclidean distance between grid locations.

Usage

imagedist(x, y, p = 2)

Arguments

Value

a list containing

distance

the Wasserstein distance $W_p(x,y)$.

plan

the optimal transport plan matrix of size $(mn\times mn)$.

Examples

#----------------------------------------------------------------------
#                       Small MNIST-like Example
#----------------------------------------------------------------------
# DATA
data(digit3)
x <- digit3[[1]]
y <- digit3[[2]]

# COMPUTE
W1 <- imagedist(x, y, p=1)
W2 <- imagedist(x, y, p=2)

# SHOW RESULTS
print(paste0("Wasserstein-1 distance: ", round(W1$distance,4)))
print(paste0("Wasserstein-2 distance: ", round(W2$distance,4)))

Description

Given two grayscale images represented as numeric matrices of identical size, compute interpolated images along a 2-Wasserstein geodesic connecting them. The function interprets each image as a discrete probability distribution on a common $(m\times n)$ grid, computes an exact optimal transport plan, and constructs intermediate measures by pushing the plan through the linear interpolation map $z=(1-t)x+t y$ (displacement interpolation / McCann's interpolation).

Usage

imageinterp(image1, image2, t = 0.5, ...)

Arguments

Details

Because the interpolated support locations generally do not coincide with the original grid points, the resulting distribution is projected back onto the grid by depositing transported mass to the nearest grid location. This is a simple and robust "re-binning" step, analogous in spirit to how histinterp re-bins interpolated quantile samples.

Value

If length(t)==1, a single $(m\times n)$ matrix representing the interpolated image. If length(t)>1, a length-length(t) list of $(m\times n)$ matrices.

Examples

#----------------------------------------------------------------------
#                  Digit Interpolation between 1 and 8
#----------------------------------------------------------------------
# LOAD DATA
set.seed(11)
data(digits)
x1 <- digits$image[[sample(which(digits$label==1),1)]] 
x2 <- digits$image[[sample(which(digits$label==8),1)]]

# COMPUTE
tvec <- seq(0, 1, length.out=10)
path <- imageinterp(x1, x2, t = tvec)

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,5), pty="s")
for (k in 1:10){
  image(path[[k]], axes=FALSE, main=sprintf("t=%.2f", tvec[k]))
} 
par(opar)

Description

Using exact balanced optimal transport as a subroutine, imagemed computes an unregularized 2-Wasserstein geometric median image $X^\dagger$ from multiple input images $X_1,\ldots,X_N$. The Wasserstein median is defined as a minimizer of the (weighted) sum of Wasserstein distances,

$\arg\min_{X} \sum_{i=1}^N w_i\, W_2(X, X_i).$

Usage

imagemed(images, weights = NULL, C = NULL, ...)

Arguments

Details

Unlike Wasserstein barycenters (which minimize squared distances), the median is a robust notion of centrality. This function solves the problem with an iterative reweighted least squares (IRLS) scheme (a Wasserstein analogue of Weiszfeld's algorithm). Each outer iteration updates weights based on current distances and then solves a weighted Wasserstein barycenter problem:

$\alpha_i^{(k)} \propto \frac{w_i}{\max(W_2(X^{(k)},X_i),\delta)}, \qquad X^{(k+1)} = \arg\min_X \sum_{i=1}^N \alpha_i^{(k)}\, W_2^2(X, X_i).$

The barycenter subproblem is solved by imagebary (mirror descent with exact OT dual subgradients). Distances $W_2$ are computed by exact EMD plans under the same squared ground cost.

Value

an $(m\times n)$ matrix of the median.

Examples

## Not run: 
#----------------------------------------------------------------------
#                             MNIST Example
#
# Use 6 images from digit '8' and 4 images from digit '1'.
# The median should look closer to the shape of '8'.
#----------------------------------------------------------------------
# DATA PREP
set.seed(11)
data(digits)
dat_8 = digits$image[sample(which(digits$label==8), 6)]
dat_1 = digits$image[sample(which(digits$label==1), 4)]
dat_all = c(dat_8, dat_1)

# COMPUTE BARYCENTER AND MEDIAN
img_bary = imagebary(dat_all, maxiter=50)
img_med  = imagemed(dat_all, maxiter=50)

# VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2), pty="s")
image(img_bary, axes=FALSE, main="Barycenter")
image(img_med,  axes=FALSE, main="Median")
par(opar)

## End(Not run)

Description

This function takes a gray-scale image represented as a matrix $X$ and converts it into a discrete measure suitable for optimal transport computations in a Lagrangian framework. Pixel intensities are normalized to sum to one, and the nonzero pixels are represented as weighted points (support and weights).

Usage

img2measure(X, threshold = TRUE)

Arguments

Value

A named list containing

support

an $(M\times 2)$ matrix of coordinates for the nonzero pixels, where each row is a point $(x,y)$.

weight

a length-$M$ vector of weights corresponding to the nonzero pixels, summing to $1$.

Examples

#-------------------------------------------------------------------
#                           Description
#
# Take a digit image and compare visualization.
#-------------------------------------------------------------------
# load the data and select the first image
data(digit3)
img_matrix = digit3[[1]]

# extract a discrete measure
img_measure = img2measure(img_matrix, threshold=TRUE)
w  <- img_measure$weight
w_norm <- w / max(w)          # now runs from 0 to 1
col_scale <- gray(1 - w_norm) # 1 = white, 0 = black

# visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2), pty="s")
image(img_matrix, xaxt="n", yaxt="n", main="Image Matrix")
plot(img_measure$support, 
     col = col_scale, xlab="", ylab="",
     pch = 19, cex = 0.5, xaxt = "n", yaxt = "n",
     main = "Extracted Discrete Measure")
par(opar)

Description

The Inexact Proximal Point Method (IPOT) offers a computationally efficient approach to approximating the Wasserstein distance between two empirical measures by iteratively solving a series of regularized optimal transport problems. This method replaces the entropic regularization used in Sinkhorn's algorithm with a proximal formulation that avoids the explicit use of entropy, thereby mitigating numerical instabilities.

Let $C := \|X_m - Y_n\|^p$ be the cost matrix, where $X_m$ and $Y_n$ are the support points of two discrete distributions $\mu$ and $\nu$, respectively. The IPOT algorithm solves a sequence of optimization problems:

$\Gamma^{(t+1)} = \arg\min_{\Gamma \in \Pi(\mu, \nu)} \langle \Gamma, C \rangle + \lambda D(\Gamma \| \Gamma^{(t)}),$