Computing the Wasserstein distance
basic_compute_distance.RmdTransport cost and Wasserstein distance
For probability measures on a metric space with finite
pth moments, the Wasserstein distance is (Villani 2003)
\[ W_p(\mu,\nu) = \left(\inf_{\gamma\in\Gamma(\mu,\nu)} \int d(x,y)^p\,d\gamma(x,y)\right)^{1/p}. \]
A coupling has the two input distributions as its marginals. For finite measures, the Kantorovich problem is a linear program over nonnegative coupling entries (Kantorovitch 1958). It permits one atom’s mass to split across several targets; a deterministic transport map need not exist for arbitrary discrete inputs.
ot_distance() distinguishes the powered transport cost
from its root. Its auto method uses monotone transport for
scalar supports and network simplex otherwise. The legacy
wasserstein() interface remains available.
A small two-dimensional example
We retain the Cassini and Smiley illustration from
mlbench, with 30 points per distribution. Standardizing
each coordinate and translating the second cloud by five defines the
geometry of this example; these transformations are part of the input
specification.
data1 <- as.matrix(scale(mlbench::mlbench.cassini(n = 30)$x))
data2 <- as.matrix(scale(mlbench::mlbench.smiley(n = 30)$x))
data2[, 1] <- data2[, 1] + 5
plot(data1, col = "blue", pch = 19, cex = .7,
xlim = range(data1[, 1], data2[, 1]),
ylim = range(data1[, 2], data2[, 2]), xlab = "x", ylab = "y")
points(data2, col = "red", pch = 19, cex = .7)
Two empirical measures after the declared scaling and translation.
Each point has equal mass within its cloud.
return_plan = TRUE retains the computed coupling for
inspection and plotting.
mu <- ot_measure(data1)
nu <- ot_measure(data2)
output <- ot_distance(mu, nu, p = 2, method = "exact", return_plan = TRUE)
output
#> Transport comparison | exact | euclidean | order 2
#> Status: SUCCESS ( optimal_pivot_condition )
#> Transport cost: 25.33738
#> Wasserstein distance: 5.033625
stopifnot(output$status == "success")
c(source_residual = max(abs(rowSums(output$plan) - mu$mass)),
target_residual = max(abs(colSums(output$plan) - nu$mass)))
#> source_residual target_residual
#> 0 0The solver reports an unregularized Wasserstein distance of 5.0336 for these generated inputs. Rows and columns of the coupling follow the original source and target atom ordering. A non-success status must not be relabeled optimal merely because its plan is feasible.
image(x = seq_len(nrow(output$plan)), y = seq_len(ncol(output$plan)),
z = output$plan, xlab = "Source index", ylab = "Target index",
main = paste("Transport plan:", toupper(output$status)))
The returned optimal coupling, in original atom order.
The same plan can be drawn as flows between locations. Line widths indicate transported mass; zero entries are omitted.
P <- output$plan
plot(data1, col = "blue", pch = 19, cex = .7,
xlim = range(data1[, 1], data2[, 1]),
ylim = range(data1[, 2], data2[, 2]), xlab = "x", ylab = "y")
points(data2, col = "red", pch = 19, cex = .7)
flows <- which(P > 0, arr.ind = TRUE)
for (k in seq_len(nrow(flows))) {
i <- flows[k, 1]; j <- flows[k, 2]
lines(c(data1[i, 1], data2[j, 1]), c(data1[i, 2], data2[j, 2]),
col = "gray50", lwd = P[i, j] / max(P))
}
Transported probability mass between the two clouds.
Supply cross-distances explicitly
A distance matrix is powered by p exactly once; an
already powered cost matrix must instead use type = "cost".
The metric flag below is justified by the Euclidean construction. For an
arbitrary dissimilarity matrix, omit that declaration and interpret the
answer as a transport cost.
D <- vapply(seq_len(nrow(data2)), function(j)
sqrt(rowSums(sweep(data1, 2, data2[j, ], "-")^2)), numeric(nrow(data1)))
crossed <- ot_distance(mu, nu, p = 2,
cost = ot_cost(D, type = "distance", metric = TRUE))
stopifnot(crossed$status == "success")
c(from_coordinates = output$wasserstein_distance,
from_cross_distances = crossed$wasserstein_distance)
#> from_coordinates from_cross_distances
#> 5.033625 5.033625For existing scripts, wassersteinD(D, p = 2) continues
to accept ground distances. The common interface adds the explicit
distinction between distance and cost matrices and a declaration of
whether a metric interpretation is intended.
Numerical scale and further workflows
A dense cost or coupling matrix requires storage proportional to the
product of source and target atom counts. Setting
return_plan = FALSE avoids retaining the plan in the
result, but does not promise a matrix-free computation. Runtime depends
on the selected solver and problem; this tutorial does not establish a
general timing or complexity bound for the implementation.
Sinkhorn regularization changes the objective and returns the transport component separately from the full regularized objective. For those contracts, input conversions, and Gaussian parameter comparisons, continue with “Distribution representations and transport objectives”. The ordinary/robust summary tutorial then demonstrates feasible sets, initialization, and failures.