T4transport (version 0.1.8)

imagebary: Barycenter of Images

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, ...)

Value

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

Arguments

images

a length-\(N\) list of same-size image matrices of size \((m\times n)\).

p

an exponent for the order of the distance (default: 2). Currently, only p=2 is supported (squared ground distance cost).

weights

a weight of each image; if NULL (default), uniform weight is set. Otherwise, it should be a length-\(N\) vector of nonnegative weights.

C

an optional \((mn\times mn)\) ground cost matrix. If NULL (default), the squared Euclidean grid cost is used. Providing C allows using alternative ground costs (e.g., geodesic distances on a manifold discretization).

...

extra parameters including

abstol

stopping criterion based on \(\ell_2\) change of iterates (default: 1e-7).

init.image

an initial barycenter image (default: arithmetic mean of normalized inputs).

maxiter

maximum number of mirror descent iterations (default: 200).

step0

initial stepsize for mirror descent (default: 0.5).

stepschedule

stepsize schedule; "sqrt" uses \(\eta_t=\text{step0}/\sqrt{t}\), and "const" uses \(\eta_t=\text{step0}\) (default: "sqrt").

eps

positivity floor for the barycenter and inputs; values are truncated below eps and renormalized (default: 1e-15). Larger values can improve robustness.

smooth

optional mixing weight toward uniform distribution after each update, used to prevent near-zero support that may cause OT infeasibility (default: 1e-12). Set to 0 to disable.

clip

\(\ell_\infty\) clipping threshold for the subgradient to stabilize exponentials in the KL update (default: 50). Set to Inf to disable.

max_backtrack

maximum number of backtracking halvings of the stepsize when an OT probe fails at the proposed update (default: 8).

print.progress

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

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.

Examples

if (FALSE) {
#----------------------------------------------------------------------
#                       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)
}