Free-Support Median by IRLS — rmedIRLS • T4transport

For a collection of empirical measures $\lbrace \mu_k\rbrace_{k=1}^K$, the free-support Wasserstein median, a minimizer to the following functional $$ \mathcal{F}(\nu) = \sum_{k=1}^K w_k \mathcal{W}_2 (\nu, \mu_k ), $$ is computed using the iteratively-reweighted least-squares (IRLS) majorization idea according to You et al. (2025) . Each surrogate receives four exact-plan relocation steps, not a guarantee of a globally solved barycenter subproblem. Accepted outer steps are checked against the full objective; local stopping does not certify global optimality.

Usage

rmedIRLS(atoms, marginals = NULL, weights = NULL, num_support = 100, ...)

Arguments

atoms

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

marginals

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

weights

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

num_support

the number of support points $M$ for the barycenter (default: 100).

...

extra parameters including

abstol: stopping criterion for iterations (default: 1e-6).
maxiter: maximum number of outer iterations (default: 10).
init.support: initial numeric support matrix; otherwise a deterministic weighted systematic selection of pooled atoms is used.
candidate.mass: fixed candidate masses; default uniform.
smoothing: nonnegative median smoothing parameter (default: zero). A positive value changes the target to the sum of square roots of squared distances plus squared smoothing.
rtol: relative support-change tolerance (default: 1e-6).
stationarity_tol: relocation-residual tolerance (default: 1e-7).
alpha: relocation step fraction in (0,1] (default: one).
inner_maxiter: exact OT optimization-pivot budget (default: 100000000); zero is unlimited.
max_backtrack: maximum objective-based step halvings (default: 12).
print.progress: show objective progress (default: FALSE).

Value

A result list retaining the legacy fields below, and adding objective, status, stopping reason, initialization, controls, and numerical diagnostics. Non-success states, including an unresolved collision, issue a warning. The feasible set has moving support locations with fixed masses. Scalar point-mass medians use a weighted-median certificate; more general collisions require a metric dominance certificate or a verified descent escape, otherwise their unresolved status remains visible.

support: an $(M \times P)$ matrix of the Wasserstein median's support points.
weight: a length-$M$ vector of median's fixed weights (uniform unless supplied).
history: accepted objective values, including initialization and genuine zeros.

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.

Examples

if (FALSE) { # \dontrun{
#-------------------------------------------------------------------
#     Free-Support Wasserstein Median of Multiple Gaussians
#
# * class 1 : samples from N((0,0),  Id)
# * class 2 : samples from N((20,0), Id)
#
#  We draw 8 empirical measures of size 50 from class 1, and 
#  2 from class 2. All measures have uniform weights.
#-------------------------------------------------------------------
## GENERATE DATA
#  8 empirical measures from class 1
input_measures = vector("list", length=10L)
for (i in 1:8){
  input_measures[[i]] = matrix(rnorm(50*2), ncol=2)
}
for (j in 9:10){
  base_draw = matrix(rnorm(50*2), ncol=2)
  base_draw[,1] = base_draw[,1] + 20
  input_measures[[j]] = base_draw
}

## COMPUTE
#  compute the Wasserstein median
run_median = rmedIRLS(input_measures, num_support = 50)
#  compute the Wasserstein barycenter
run_bary   = rbaryGD(input_measures, num_support = 50)

## VISUALIZE
opar <- par(no.readonly=TRUE)

#  draw the base points of two classes
base_1 = matrix(rnorm(80*2), ncol=2)
base_2 = matrix(rnorm(20*2), ncol=2)
base_2[,1] = base_2[,1] + 20
base_mat = rbind(base_1, base_2)
plot(base_mat, col="gray80", pch=19)

#  auxiliary information
title("estimated barycenter and median")
abline(v=0); abline(h=0)

#  draw the barycenter and the median
points(run_bary$support, col="red", pch=19)
points(run_median$support, col="blue", pch=19)
par(opar)
} # }