Ordinary and robust summaries with explicit numerical status

Compare targets on the same observations

A barycenter minimizes a weighted sum of squared W2 distances. A median minimizes a weighted sum of W2 distances, so its numerical objective is in different units. Compare their fitted distributions and per-observation residuals; a smaller raw objective across these different targets is not a model-selection criterion.

observations <- ot_collection(lapply(c(0, 0, 0, 2, 8), ot_measure))
ordinary <- ot_summary(observations, target = "barycenter")
robust <- ot_summary(observations, target = "median")
c(barycenter_location = ordinary$estimate$support[1, 1],
  median_location = robust$estimate$support[1, 1])
#> barycenter_location     median_location 
#>                   2                   0
summary(robust)$per_measure
#>      label weight distance contribution
#> 1 measure1    0.2        0          0.0
#> 2 measure2    0.2        0          0.0
#> 3 measure3    0.2        0          0.0
#> 4 measure4    0.2        2          0.4
#> 5 measure5    0.2        8          1.6

In this deliberately simple example, the majority of input distributions are point masses at zero. The barycenter lies at two; the median stays at zero. This analytical illustration is not evidence of robustness for every data representation, contamination pattern, or local solver.

State what is allowed to change

Feasible set	Optimized quantities	Main limitation
unrestricted_1d	Full scalar quantile function on its common refinement	Scalar measures only
free_support_fixed_mass	A fixed number of candidate locations	Candidate masses stay fixed; optimization is generally local
fixed_support	Probability masses on supplied locations	Estimate cannot leave the supplied grid

Here is a small two-dimensional collection of translated four-atom measures. Three of the five observations coincide. We provide the same support and candidate masses to both fits so the comparison has the same feasible set.

shape <- rbind(c(0, 0), c(1, 0), c(0, 1), c(1, 1))
shifts <- c(0, 0, 0, 2, 8)
clouds <- lapply(shifts, function(s) ot_measure(sweep(shape, 2, c(s, 0), "+")))
ordinary_cloud <- ot_summary(clouds, target = "barycenter",
  feasible = "free_support_fixed_mass", support = shape, mass = rep(1/4, 4))
robust_cloud <- ot_summary(clouds, target = "median",
  feasible = "free_support_fixed_mass", support = shape, mass = rep(1/4, 4))
stopifnot(ordinary_cloud$status == "success", robust_cloud$status == "success")
rbind(barycenter_center = colSums(ordinary_cloud$estimate$support * ordinary_cloud$estimate$mass),
      median_center = colSums(robust_cloud$estimate$support * robust_cloud$estimate$mass))
#>                   [,1] [,2]
#> barycenter_center  2.5  0.5
#> median_center      0.5  0.5
robust_cloud$diagnostics$optimality_certificate
#> [1] "metric median with at least half the outer mass at the returned measure"

The barycenter is the shape translated by the average shift; the median is the majority shape. The reported median certificate concerns this feasible majority input. More general fits need not have such a certificate.

op <- par(mfrow = c(1, 2))
plot(ordinary_cloud$estimate, main = "Barycenter", xlim = c(-.5, 3.5), ylim = c(-.5, 1.5))
plot(robust_cloud$estimate, main = "Median", xlim = c(-.5, 3.5), ylim = c(-.5, 1.5))

The ordinary summary follows the mean translation; the median retains the majority shape.

par(op)

Initialization and support-size sensitivity

Without supplied locations, the first free-support start uses deterministic systematic selection from pooled input atom indices weighted by both outer and within-measure masses. Later starts sample pooled atoms with those probabilities. This initialization depends on atom order. Supply support when the initial geometry must be controlled.

multiple <- ot_summary(clouds, target = "barycenter",
  feasible = "free_support_fixed_mass", support = shape, mass = rep(1/4, 4),
  nstart = 3, seed = 23)
data.frame(start = seq_len(multiple$multistart$nstart),
           objective = multiple$multistart$objectives,
           status = multiple$multistart$statuses,
           reason = multiple$multistart$stop_reasons)
#>   start objective  status              reason
#> 1     1       9.6 success relocation_residual
#> 2     2       9.6 success relocation_residual
#> 3     3       9.6 success relocation_residual
multiple$multistart$chosen
#> [1] 1
multiple$initialization
#> $support
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    1    0
#> [3,]    0    1
#> [4,]    1    1
#> 
#> $mass
#> [1] 0.25 0.25 0.25 0.25
#> 
#> $method
#> [1] "supplied support"

The selected run has the smallest objective among successful runs. If every run fails, the best finite failed result is retained with a non-success status. The local seed restores the caller’s prior RNG state. Saving only a seed is insufficient for an audit: retain the initial supports, masses, controls, statuses, and fitted objective as well.

Increasing the candidate support size changes the feasible approximation. The table below reports both the objective and numerical outcome; a failed run must not be silently dropped from a sensitivity analysis.

sizes <- c(1L, 2L, 4L)
size_fits <- lapply(sizes, function(m) ot_summary(clouds, target = "barycenter",
  feasible = "free_support_fixed_mass", control = list(num_support = m),
  nstart = 3, seed = 23))
data.frame(atoms = sizes,
           objective = vapply(size_fits, function(x) x$objective, numeric(1)),
           status = vapply(size_fits, function(x) x$status, character(1)))
#>   atoms objective  status
#> 1     1     10.10 success
#> 2     2      9.85 success
#> 3     4      9.60 success

Read failures and stopping criteria

history contains accepted objective values, starting with initialization. The plot and printed result show status. Relocation residuals use the currently selected optimal plans; a small residual does not generally certify global optimality over moving supports.

limited <- ot_summary(clouds, target = "barycenter",
  feasible = "free_support_fixed_mass", support = shape, mass = rep(1/4, 4),
  control = list(maxiter = 1, alpha = 0.1))
#> Warning: Summary computation: iteration_limit (maximum_iterations).
limited[c("status", "stop_reason", "objective")]
#> $status
#> [1] "iteration_limit"
#> 
#> $stop_reason
#> [1] "maximum_iterations"
#> 
#> $objective
#> [1] 12.84
plot(limited, type = "history")

The warning above is intentional: a single damped update does not meet the relocation criterion. An inner_solver_failure instead means a required transport subproblem did not succeed. A step_rejected outcome means no acceptable trial update was established. An unresolved collision is a failure condition, not proof that the current observation is a median. Increasing an iteration budget may address a limit; it is not a universal repair for other failure reasons.

A positive median control$smoothing changes the target to sum(weights * sqrt(W2^2 + smoothing^2)), reported as smoothed_median. This is different from entropy regularization of an inner transport plan. Compare runs only when target, smoothing, candidate masses, support size, and regularization are held fixed or explicitly part of the sensitivity study.

Fixed-grid summaries

On a prescribed grid, the candidate masses vary instead of the locations. Initial masses use control$init_mass and must be strictly positive for the mirror-descent implementation. Input measures themselves may contain zeros.

fixed <- ot_summary(list(ot_measure(0), ot_measure(3)), weights = c(.75, .25),
  target = "barycenter", feasible = "fixed_support", support = 0:3)
fixed
#> Wasserstein barycenter | fixed_support | subgradient 
#> Status: SUCCESS ( convex_dual_bound_tolerance )
#> Objective: 1.75 | estimated atoms: 4
fixed$diagnostics[c("convex_lower_bound", "objective_gap")]
#> $convex_lower_bound
#> [1] 1.75
#> 
#> $objective_gap
#> [1] 2.606099e-08
regularized <- ot_summary(list(ot_measure(0)), target = "barycenter",
  feasible = "fixed_support", support = 0:1, epsilon = 0.5)
regularized
#> Wasserstein barycenter | fixed_support | subgradient 
#> Status: SUCCESS ( single_target_columnwise_entropy_optimum )
#> Objective: -0.563464 | estimated atoms: 2 
#> Entropy coefficient: 0.5 (full regularized objective)

The unregularized barycenter objective is convex in the candidate mass vector; its reported valid dual lower bound can support an objective-gap stopping rule. The fixed-grid median objective is generally nonconvex and has no corresponding global convex-gap claim. Regularized median requests are rejected. The regularized barycenter above optimizes the complete entropy objective, including the candidate-mass entropy contribution; it is not the unregularized solution.

These tutorials run offline and use only small deterministic examples. Full application studies, data provenance records, and discretization sensitivity results are maintained in the manuscript’s separate replication bundle.