Package 'T4transport'

Description

digit3 contains 2000 images from the MNIST dataset of digit 3. Each element of the list is an image represented as an $(28\times 28)$ matrix that sums to 1. Normalization makes intensities into probability masses and discards the original total intensity of each image.

Usage

data(digit3)

Format

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

Source

MNIST Database of Handwritten Digits, by Yann LeCun, Corinna Cortes, and Christopher J.C. Burges, http://yann.lecun.com/exdb/mnist/. The original package source does not preserve this subset's MNIST row indices, train/test designation, random seed, or exact preparation code. See the installed ‘DATA_PROVENANCE’ for source-artifact checksums and the limits of the provenance and dataset-specific licensing evidence.

References

LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998). Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11), 2278–2324. doi:10.1109/5.726791.

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 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. Normalization makes intensities into probability masses and discards the original total intensity of each image.

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.

Source

MNIST Database of Handwritten Digits, by Yann LeCun, Corinna Cortes, and Christopher J.C. Burges, http://yann.lecun.com/exdb/mnist/. The retained preparation fragment uses dslabs::read_mnist(), samples 500 training images per class without replacement, reshapes each row with matrix(pixel_vector, nrow = 28)[, 28:1], and normalizes total mass. Its original random seed, selected MNIST row indices, and dslabs version were not recorded. See the installed ‘DATA_PROVENANCE’ for checksums, transformations, and limits of the dataset-specific licensing evidence.

References

LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998). Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11), 2278–2324. doi:10.1109/5.726791.

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 on their exact common refinement. Ties and unequal empirical sample sizes retain their probability masses.

Usage

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

Arguments

Value

A callable, right-continuous weighted CDF of class "t4transport_cdf", which also inherits from "stepfun". Its "t4transport_measure" attribute contains exact support locations and masses; "t4transport_fit" contains the objective, accepted history, and numerical diagnostics. It can be supplied to subsequent calls to ecdfbary() or ecdfmed(). No artificial repeated sample is constructed.

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 collision-aware modified Weiszfeld algorithm in exact finite quantile coordinates. The target is a sum of Wasserstein-2 distances, not coordinatewise quantile medians. A warning is issued if the stationarity criterion is not met; details remain in the returned fit attribute.

Usage

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

Arguments

Value

A callable, right-continuous weighted CDF of class "t4transport_cdf", which also inherits from "stepfun". Its "t4transport_measure" attribute contains exact support locations and masses; "t4transport_fit" contains the objective, accepted history, and numerical diagnostics. It can be supplied to subsequent calls to ecdfbary() or ecdfmed(). No artificial repeated sample is constructed.

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. The compatibility implementation is described in Details.

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

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.

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.

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 = 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 = 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. 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

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.

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.

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

Details

This graph summary is not a transport distance. The eigenvector sign is fixed by making its largest absolute component positive; repeated eigenvalues can still give different bases on different platforms. At least two points and a positive median pairwise distance are required.

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-2 barycenter exactly: average the means and standard deviations separately, then square the resulting standard deviation. The iteration controls are retained and validated for compatibility.

Usage

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

Arguments

Value

A named list preserving the mean and variance fields below, with objective, accepted history, status, stopping reason, iteration count, and numerical diagnostics. Non-success produces a warning.

mean

mean of the estimated barycenter distribution.

var

variance of the estimated barycenter distribution.

References

There are no references for Rd macro ⁠\insertAllCites⁠ on this help page.

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 $p$-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 preserving the mean and variance fields below, with objective, accepted history, status, stopping reason, iteration count, and numerical diagnostics. Non-success produces a warning.

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 median under Wasserstein-2 distance within the Gaussian family. This is a collision-aware Euclidean geometric median of the pairs (mean, standard deviation), not separate coordinatewise medians.

Usage

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

Arguments

Value

A named list preserving the mean and variance fields below, with objective, accepted history, status, stopping reason, iteration count, and numerical diagnostics. Non-success produces a warning.

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-2 median within the strictly SPD Gaussian family. General covariance fits use iteratively reweighted squared-distance subproblems with arithmetic mean updates and covariance barycenters. Objective and iterate tolerances are numerical stopping criteria; they do not certify a globally optimal median.

Usage

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

Arguments

Value

A named list preserving the mean and variance fields below, with objective, accepted history, status, stopping reason, iteration count, and numerical diagnostics. Non-success produces a warning.

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 places deterministic 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

Details

The contour has Mahalanobis radius one; it is not a random sample or an automatically calibrated confidence region. Covariance must be strictly positive definite, with no silent ridge or projection repair.

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 symmetric dissimilarity space representing the input collection under a local GW squared-loss criterion. Triangle inequalities are not imposed, so the returned object need not be a metric space. Global optimality is not certified.

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 attained dissimilarity barycenter; a metric is not guaranteed.

weight

a length-$M$ vector of fixed barycenter masses with all entries being $1/M$. Status, objective history, dissimilarity matrix, and pairwise diagnostics are also returned.

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(digit3)
data_D = vector("list", length=5)
data_W = vector("list", length=5)
for (i in 1:5){
  img_now = img2measure(digit3[[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 a local Gromov-Wasserstein (GW) comparison between finite spaces. Given two distance matrices $D_X$ and $D_Y$ along with their respective marginal distributions, the function solves the GW optimization problem approximately to obtain a squared-loss objective and an associated coupling. Global optimality is not certified.

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: Frank–Wolfe (FW) with quadratic line search, a full linear-oracle update (the historical MM label), and projected gradient (PG) with Dykstra Euclidean projection and backtracking. These are unregularized nonconvex iterations. Symmetry and zero-diagonal checks do not establish the triangle inequality. The objective is the full squared-loss sum, without a factor 1/2.

Usage

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

Arguments

Details

Squared-loss evaluations use a fast matrix expansion, with direct nonnegative block sums where a roundoff bound indicates cancellation. This avoids silently reporting zero for nearly equal large dissimilarities. The fallback uses bounded working memory but may require quartic arithmetic in the worst case. Costs outside floating-point range fail explicitly.

Value

a named list containing

distance

square root of the attained squared-loss objective, not a certified global GW distance.

plan

an $(M\times N)$ nonnegative coupling matrix. The result also records objective, status, first-order gap, and accepted history.

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

Compute the finite 1D 2-Wasserstein barycenter of midpoint-discrete histograms, then return its projection onto the original histogram bins. All input histograms must have identical breaks. The unprojected finite measure is retained separately.

Usage

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

Arguments

Details

Inputs are discrete measures at bin midpoints, using normalized counts (or the canonical mass field of an estimated histogram). The unprojected finite result is saved in unprojected_measure. For the legacy histogram display it is binned onto the original breaks, using right-closed bins with the lowest endpoint included, and represented at bin midpoints. This projection generally changes the measure; it is not an exact unrestricted Wasserstein summary or geodesic. Downstream histogram functions consistently use the displayed midpoint measure, not the separate unprojected result. Break values must agree exactly across inputs; no grid snapping is done. When an estimated histogram is reused, its previous computational metadata is discarded. Returned diagnostics describe the current calculation. Returned mass and counts are the same probabilities; density equals mass divided by bin width. Counts are not rounded observations, and effective_sample_size is NA because an estimated measure has no implied sample size. The objective/history and solver status refer to the unprojected summary.

Value

a "histogram" display of the barycenter, with the unprojected estimate and projection metadata.

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. Break values must agree exactly across inputs; no grid snapping is done.

Usage

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

Arguments

Details

Each input is the discrete measure at its bin midpoints with normalized bin counts, including unequal-width bins. The entropy coefficient lambda multiplies the sum of transport-plan entries times their logarithms; the cost is midpoint distance raised to p. The finite regularization coefficient changes the objective from the unregularized barycenter target. The shared solver now uses log-domain Sinkhorn gradients and objective-checked mirror descent; historical names remain compatibility entry points. Returned mass and counts are probability masses, and density is mass divided by bin width. No sample-size-dependent rounding is performed; effective_sample_size is NA. Status, accepted objective history, and inner diagnostics are available in fit. The full objective uses log(plan) - 1; this differs from log(plan) by a constant and has the same minimizer.

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. Break values must agree exactly across inputs; no grid snapping is done.

Usage

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

Arguments

Details

Each input is the discrete measure at its bin midpoints with normalized bin counts, including unequal-width bins. The entropy coefficient lambda multiplies the sum of transport-plan entries times their logarithms; the cost is midpoint distance raised to p. The finite regularization coefficient changes the objective from the unregularized barycenter target. The shared solver now uses log-domain Sinkhorn gradients and objective-checked mirror descent; historical names remain compatibility entry points. Returned mass and counts are probability masses, and density is mass divided by bin width. No sample-size-dependent rounding is performed; effective_sample_size is NA. Status, accepted objective history, and inner diagnostics are available in fit. The full objective uses log(plan) - 1; this differs from log(plan) by a constant and has the same minimizer.

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 or canonical estimated bin masses. Uses the exact 1D monotone OT algorithm, not LP nor entropic regularization. Break values must agree exactly; to compare different grids, convert each histogram with as_ot_measure() and use ot_distance().

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

Compute exact finite displacement interpolation between midpoint-discrete histograms with identical breaks, then project each result onto the original histogram bins. The unprojected finite measures follow the 2-Wasserstein geodesic; the projected histogram displays generally do not.

Usage

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

Arguments

Details

Inputs are discrete measures at bin midpoints, using normalized counts (or the canonical mass field of an estimated histogram). The unprojected finite result is saved in unprojected_measure. For the legacy histogram display it is binned onto the original breaks, using right-closed bins with the lowest endpoint included, and represented at bin midpoints. This projection generally changes the measure; it is not an exact unrestricted Wasserstein summary or geodesic. Downstream histogram functions consistently use the displayed midpoint measure, not the separate unprojected result. Break values must agree exactly across inputs; no grid snapping is done. When an estimated histogram is reused, its previous computational metadata is discarded. Returned diagnostics describe the current calculation. Returned mass and counts are the same probabilities; density equals mass divided by bin width. Counts are not rounded observations, and effective_sample_size is NA because an estimated measure has no implied sample size. The objective/history and solver status refer to the unprojected summary.

Value

If length(t) == 1, a single "histogram" object representing the projected interpolation at time t, including the unprojected measure. 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

Compute a 1D Wasserstein geometric median of midpoint-discrete histograms with common breaks, using a modified Weiszfeld algorithm on exact finite quantile intervals. Return its projection onto the original histogram bins and retain the unprojected finite estimate separately.

Usage

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

Arguments

Details

Inputs are discrete measures at bin midpoints, using normalized counts (or the canonical mass field of an estimated histogram). The unprojected finite result is saved in unprojected_measure. For the legacy histogram display it is binned onto the original breaks, using right-closed bins with the lowest endpoint included, and represented at bin midpoints. This projection generally changes the measure; it is not an exact unrestricted Wasserstein summary or geodesic. Downstream histogram functions consistently use the displayed midpoint measure, not the separate unprojected result. Break values must agree exactly across inputs; no grid snapping is done. When an estimated histogram is reused, its previous computational metadata is discarded. Returned diagnostics describe the current calculation. Returned mass and counts are the same probabilities; density equals mass divided by bin width. Counts are not rounded observations, and effective_sample_size is NA because an estimated measure has no implied sample size. The objective/history and solver status refer to the unprojected summary.

Value

a "histogram" display of the median, with the unprojected estimate, solver diagnostics, and projection metadata.

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 optional candidate smoothing and backtracking against the full objective. Input masses are never floored. Inner-solver failures return a failed result; they do not trigger an automatic change to the input distributions.

Value

An $(m\times n)$ matrix with a diagnostics attribute, or a diagnostic result list when return_result = TRUE. Non-successful matrix returns issue a warning. Small iterate change does not certify global optimality of this nonsmooth numerical optimization.

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

Details

Since version 0.2.0 the historical names call the shared objective-checked mirror-descent solver with log-domain Sinkhorn gradients, as documented in fbary14C. The entropy coefficient lambda multiplies sum(plan * (log(plan) - 1)). Pixel columns run from x=0 to x=1 and rows from y=1 to y=0; a singleton axis is fixed at its first coordinate. Pixel intensities are probability masses after normalization. Additional controls from fbary14C are accepted, using init.image instead of init.vec.

Value

A matrix of image probability masses with a diagnostics attribute; return_result = TRUE returns a rich fit and 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

Details

Since version 0.2.0 the historical names call the shared objective-checked mirror-descent solver with log-domain Sinkhorn gradients, as documented in fbary14C. The entropy coefficient lambda multiplies sum(plan * (log(plan) - 1)). Pixel columns run from x=0 to x=1 and rows from y=1 to y=0; a singleton axis is fixed at its first coordinate. Pixel intensities are probability masses after normalization. Additional controls from fbary14C are accepted, using init.image instead of init.vec.

Value

A matrix of image probability masses with a diagnostics attribute; return_result = TRUE returns a rich fit and 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
datsmall = list(matrix(c(1, 0, 0, 1), 2, 2))
outsmall = imagebary15B(datsmall, lambda=0.5, maxiter=3)
attr(outsmall, 'diagnostics')$status

## 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, then projects it onto the original grid. 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. The unprojected finite measure is retained in an attribute; the displayed matrix generally is not an exact Wasserstein geodesic.

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

Minimizes the weighted sum of W2 distances over probability masses on the common image grid, using exact OT values and objective-checked mass updates. The mass-space median objective is nonconvex. A local stopping criterion is not a global-optimality certificate; a dominating input can be certified separately by the triangle inequality for the default Euclidean geometry.

Usage

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

Arguments

Details

Since version 0.2.0, the old denominator-floor IRLS and its nested bary.*/delta controls have been replaced by the shared fixed-support median solver. Input zero masses are retained exactly. An uncertified stationary point, ambiguous zero-distance collisions, and inner solver failures are reported; no failed input is discarded or downweighted. Pixel columns span [0,1] and rows run from y=1 to y=0, as in imagedist.

Value

An image probability-mass matrix with a diagnostics attribute, or a full fit when return_result = TRUE.

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

Details

This legacy converter uses pixel units and base image()'s row-as-x orientation. It differs from the endpoint grid of imagedist and the pixel-centre grid of ot_image. Choose a coordinate convention explicitly when comparing measures. Thresholding reports removed mass.

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 KL proximal formulation. Log scaling reduces kernel underflow but does not remove all floating-point limitations.

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)}),$

where $\lambda > 0$ is the proximal regularization parameter and $D(\cdot \| \cdot)$ is the Kullback–Leibler divergence. Each subproblem is solved approximately using a fixed number of inner iterations, making the method inexact.

The returned value is the rooted transport cost of an approximate plan. Small changes in the plan and small marginal residuals are stopping criteria, not a certificate of optimality for the unregularized problem. The parameter lambda controls proximal steps and is not an entropy penalty in the reported objective.

Usage

ipot(X, Y, p = 2, wx = NULL, wy = NULL, lambda = 1, ...)

ipotD(D, p = 2, wx = NULL, wy = NULL, lambda = 1, ...)

Arguments