Distance-Matrix Wasserstein Statistics for Scalable Gromov–Wasserstein Learning
arXiv:2605.14981v1 [cs.LG] 14 May 2026
Ao Xu
School of Artificial Intelligence Jilin University No. 2699, Qianjin Street, Chaoyang District Changchun 130012, China and Zhongguancun Academy Daniufang 2nd Ring Road, Haidian District Beijing 100094, China
Tieru Wu∗
School of Artificial Intelligence Jilin University No. 2699, Qianjin Street, Chaoyang District Changchun 130012, China
Abstract Gromov–Wasserstein (GW) distances compare graphs, shapes, and point clouds through internal distances, without requiring a common coordinate system. This invariance is powerful, but discrete GW is a nonconvex quadratic optimal transport problem and is difficult to estimate at scale. We propose Distance-Matrix Wasserstein (DMW), a hierarchy of Wasserstein statistics comparing laws of random finite distance matrices. Rather than optimizing a global point-level alignment, DMW samples n points from each space, records their pairwise distances, and transports the resulting matrix laws. We prove that DMW is a relaxation and lower bound of GW, and establish a reverse approximation inequality: the GW–DMW gap is controlled by the Wasserstein error of approximating each original measure with n samples. Hence population DMW converges to GW as sampled subspaces become dense. We further give finite-sample bounds, including intrinsic-dimensional rates that depend on the data manifold rather than the ambient matrix dimension n2 . For scalable computation, we introduce sliced and multi-scale DMW; for p = 1, the sliced multiscale dissimilarity yields positive-definite exponential kernels. Experiments on synthetic metric spaces, scalability benchmarks, graph classification, and two-sample testing validate the theory and demonstrate an interpretable GW-style proxy for structural comparison. Keywords: Gromov–Wasserstein distance, metric measure spaces, distance matrix distributions, sliced optimal transport, graph kernels, scalable structured learning
1 Introduction Many modern learning problems compare objects whose points do not live in a shared coordinate system. A molecule, a protein interaction network, a social graph, and a point cloud sampled from a shape can all be viewed as structured objects, but their vertex sets, labels, and ambient spaces are usually unrelated. In such settings, similarity is not primarily ∗. Corresponding author. ©2026 Ao Xu and Tieru Wu. License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/.
Xu and Wu
a question of Euclidean position. It is a question of internal geometry: do neighborhoods, paths, clusters, and long-range distances inside one object resemble those inside another? Gromov–Wasserstein (GW) distances give one of the most principled answers to this question (Gromov, 1999; Mémoli, 2011; Villani, 2009; Peyré and Cuturi, 2019). Rather than comparing coordinates, GW compares pairwise distance patterns. It looks for a coupling between the points of two metric measure spaces such that pairs that are close in one space are matched, on average, to pairs that are close in the other. This intrinsic viewpoint has made GW and its variants useful for graphs, shapes, point clouds, and structured optimal transport problems (Vayer et al., 2019a). The same feature, however, is also the main bottleneck. Discrete GW optimizes over a global point-level coupling and has a nonconvex quadratic objective. For large graphs or point clouds, this global alignment is expensive to compute, unstable to estimate, and often more detailed than what is needed for downstream tasks such as classification, retrieval, or two-sample testing. There is a long history of replacing hard structural alignment problems by more tractable summaries. Graph kernels compare shortest-path statistics, graphlets, or Weisfeiler–Lehman refinement features (Borgwardt and Kriegel, 2005; Shervashidze et al., 2009, 2011). Distance distributions and related shape statistics provide lower-dimensional invariants for metric spaces and are closely connected to lower bounds for GW (Mémoli, 2011; Mémoli and Needham, 2022). In a different direction, sliced and projected optimal transport methods reduce high-dimensional transport to many low-dimensional problems (Bonneel et al., 2015; Vayer et al., 2019b; Lin et al., 2020). These ideas are computationally attractive, but they also raise a statistical question: can one design a scalable statistic that is still tied quantitatively to GW, rather than merely being a convenient feature? The starting point of this paper is the observation that much of a metric space is visible through random finite snapshots. Imagine sampling n points from a graph, a point cloud, or a shape, and recording all pairwise distances among the sampled points. The resulting distance matrix is a small finite metric subspace, a local-to-global fingerprint of the object. Repeating this experiment produces a probability law on distance matrices. If two objects generate similar laws of such finite metric snapshots, then their random finite substructures look similar even without constructing a global pointwise alignment. This viewpoint is classical at the level of invariants. The law of the infinite random distance array determines a metric measure space up to measure-preserving isometry, a principle going back to Gromov and Vershik and developed probabilistically through the Gromov-weak topology (Gromov, 1999; Vershik, 2004; Greven et al., 2009). Our goal is different. We do not use distance matrix distributions only as reconstruction objects. We use their finite-dimensional laws as computable statistics, and we ask how their Wasserstein comparison relates to GW. We call the resulting statistic Distance-Matrix Wasserstein (DMW). At order n, DMW samples n points from each space, maps each sample to the vector of its pairwise distances, and computes a Wasserstein distance between the two induced laws on distance matrices. Conceptually, DMW relaxes GW by moving the coupling from the point level to the finitesubspace level. GW requires a single coupling that consistently aligns all sampled pairs through the same underlying point matching. DMW instead transports random distance matrices directly. This removes the explicit alignment, but it keeps the comparison intrinsic and makes it naturally sample based. Figure 1 illustrates the computational distinction: 2
Distance-Matrix Wasserstein Statistics
GW
DMW
global point-level alignment X
transport distance-matrix laws
Y X
DX DX DX
ρX
Y
DY DY DY
(n) ρY OT on matrix laws
sample n-tuples
one coupling π all pairwise distortions are evaluated through the same point matching
distance matrices
(n)
Wp transport laws
Figure 1: The computational contrast between GW and DMW. GW optimizes a global point-level coupling and evaluates pairwise distance distortion through that coupling. DMW samples finite metric subspaces, maps each sample to a distance matrix, and performs optimal transport between the resulting distance-matrix laws.
GW couples original points and evaluates pairwise distortions through one global matching, whereas DMW first converts random finite subspaces into distance matrices and then transports the induced matrix laws. The main theoretical message is that this relaxation is controlled. We prove first that DMW is always a lower bound for GW. More importantly, we prove a reverse approximation inequality: the gap between GW and order-n DMW is bounded by the error of approximating each original measure by an empirical measure with n samples. In plain terms, DMW is close to GW whenever n random points are enough to represent the two spaces in Wasserstein distance. Thus, at the population level, increasing the order makes the finitesubspace statistic converge back to GW under standard empirical Wasserstein convergence assumptions (Fournier and Guillin, 2015). This result is the central bridge of the paper: DMW is not only a lower-bound feature, but a consistent finite-subspace approximation to GW. The same theorem also exposes the central tradeoff. Larger orders contain richer geometry and improve the population approximation, but the distance-matrix law lives in dimension n(n − 1)/2 and becomes harder to estimate directly. We therefore develop empirical, sliced, and multi-scale versions of DMW. Slicing replaces high-dimensional optimal transport between distance-matrix laws by one-dimensional Wasserstein computations along random directions, while the multi-scale statistic averages several orders to balance local distance information, higher-order geometry, and finite-sample stability. The theory makes this tradeoff explicit through concentration bounds, finite-direction estimates, intrinsicdimensional rates, and a positive-definiteness result for the p = 1 sliced kernel used in graph experiments. The empirical study is designed to test this story rather than only to report benchmark accuracy. Synthetic experiments check the approximation–estimation tradeoff predicted by the theory and the effect of the number of slicing directions. Scalability experiments isolate 3
Xu and Wu
the dependence on the order n, the number of sampled matrices K, and the number of directions L. Graph classification and two-sample testing experiments place the method in the regime where it is intended to be useful: unsupervised, intrinsic, GW-style structural comparison without solving a global alignment problem for every pair of large objects. Contributions.
We summarize the contributions in four parts.
1. We introduce DMW, a hierarchy of Wasserstein distances between random finite distancematrix laws, and show that it is an intrinsic relaxation and lower bound of GW. 2. We establish the main reverse approximation theorem: the GW–DMW gap is controlled by empirical Wasserstein approximation errors. This yields population consistency, finite-space rates, Euclidean rates, and manifold/intrinsic-dimensional samplecomplexity bounds. 3. We develop empirical, sliced, and multi-scale DMW estimators, including high-probability concentration, finite-direction slicing guarantees, sliced expressivity via Cramer–Wold arguments, and a positive-definite p = 1 kernel for downstream learning. 4. We provide a reproducible experimental suite validating the theory on synthetic metric spaces and evaluating scalability, two-sample testing, and graph classification on ten TU benchmarks.
2 Preliminaries This section fixes the notation and conventions used throughout the paper. We work with compact metric measure spaces in the main text in order to avoid moment and measurability distractions. The arguments extend to complete separable metric spaces with finite pmoments whenever the relevant disintegrations and Wasserstein distances are well-defined. 2.1 Notation We use calligraphic letters such as X = (X, dX , µ) and Y = (Y, dY , ν) for metric measure spaces, and plain letters X, Y for their underlying compact metric spaces. The symbols µ, ν always denote the corresponding probability measures. For a measurable map T , the pushforward of a measure α is written T# α, and Π(α, β) denotes the set of couplings of α and β. The ordinary Wasserstein distance on a metric space Z is denoted by Wp,Z , or simply Wp when the underlying space is clear. The Gromov–Wasserstein distance between metric measure spaces is denoted by GWp . For an integer n ≥ 2, we write n Nn = 2 (n)
for the number of unordered pairs in an n-tuple. The map TX : X n → RNn sends a tuple to its vector of pairwise distances, and (n)
(n)
ρX = (TX )# µ⊗n 4
Distance-Matrix Wasserstein Statistics
is the corresponding distance-matrix law. The population, empirical, sliced, multi-scale sliced, and unsliced multi-scale statistics are denoted respectively by DMWn,p ,
\ n,p,K , DMW
SDMWn,p ,
MSDMWp ,
MDMWp .
The parameter n is the finite-subspace order, K is the number of sampled distance matrices, and L is the number of slicing directions. Covering numbers are written as N (A, d, ε). When a dimension parameter is needed for a covering bound, we use sX , sY rather than dX , dY , reserving dX , dY for metrics. 2.2 Probability Measures and Couplings Let (Z, dZ ) be a compact metric space and let B(Z) denote its Borel σ-algebra. We write P(Z) for the set of Borel probability measures on Z. For α, β ∈ P(Z), a coupling of α and β is a probability measure γ ∈ P(Z × Z) satisfying (pr1 )# γ = α,
(pr2 )# γ = β,
where pr1 , pr2 are the coordinate projections. The set of all such couplings is denoted by Π(α, β). More generally, if α ∈ P(Z) and β ∈ P(Z ′ ) are defined on possibly different compact metric spaces, then Π(α, β) denotes the set of probability measures on Z ×Z ′ with marginals α and β. If T : Z → Z ′ is measurable, T# α denotes the pushforward measure. 2.3 Wasserstein Distance For p ≥ 1, the p-Wasserstein distance between α, β ∈ P(Z) is Wp,Z (α, β) =
Z
1/p dZ (z, z ) dγ(z, z ) . ′ p
inf γ∈Π(α,β) Z×Z
′
When the underlying metric space is clear, we write Wp instead of Wp,Z . Compactness ensures that the infimum is finite and attained. We will repeatedly use the elementary inequality |dZ (z, z ′ ) − dZ (w, w′ )| ≤ dZ (z, w) + dZ (z ′ , w′ ), which follows from the triangle inequality. 2.4 Metric Measure Spaces A metric measure space is a triple X = (X, dX , µ), where (X, dX ) is a compact metric space and µ ∈ P(X). Two metric measure spaces X = (X, dX , µ) and Y = (Y, dY , ν) are measure-preserving isometric if there exists an isometry φ : supp(µ) → supp(ν) such that φ# µ = ν. Throughout the paper, statements about equality of metric measure spaces are understood up to this equivalence. 5
Xu and Wu
Given independent samples x1 , . . . , xn ∼ µ, we write n
1X µ bn = δxi , n
Xbn = (X, dX , µ bn ).
i=1
The quantity EWp,X (µ, µ bn ) measures how well n random samples approximate µ in the intrinsic metric of X. 2.5 Gromov–Wasserstein Distance For compact metric measure spaces X = (X, dX , µ) and Y = (Y, dY , ν), we use the convention ZZ p GWp (X , Y) = inf |dX (x, x′ ) − dY (y, y ′ )|p dπ(x, y) dπ(x′ , y ′ ). π∈Π(µ,ν)
(X×Y )2
The coupling π aligns mass between the two spaces, and the objective measures the average distortion of pairwise distances under this alignment. With this convention, GWp is invariant under measure-preserving isometries and descends to a metric on the corresponding quotient class of compact metric measure spaces. Remark 2.1 (Normalization convention) Some authors include constants such as 1/2 or 1/21/p in the definition of GWp . We do not include such constants. Consequently, stability constants below are calibrated for the definition above. 2.6 Finite Distance-Matrix Spaces For an integer n ≥ 2, let Nn =
n n(n − 1) = 2 2
and identify RNn with arrays indexed by unordered pairs (i, j) with 1 ≤ i < j ≤ n. For a, b ∈ RNn , define the normalized ℓp metric 1/p X 1 ∥a − b∥p,avg = |aij − bij |p . Nn 1≤i<j≤n
The normalization makes the scale of the distance independent of the number of pairwise entries. For a metric measure space X = (X, dX , µ), define the order-n distance-matrix map (n)
TX : X n → RNn , (n)
Since dX is continuous, TX matrix law is
(n)
TX (x1 , . . . , xn ) = (dX (xi , xj ))1≤i<j≤n .
is Borel measurable. The corresponding random distance(n)
(n)
ρX = (TX )# µ⊗n ∈ P(RNn ). (n)
Equivalently, ρX is the law of the pairwise distance vector generated by n independent samples from µ. 6
Distance-Matrix Wasserstein Statistics
3 Distance-Matrix Wasserstein This section introduces the object that drives the paper. The guiding idea is simple: rather than compare two spaces by directly aligning their points, compare the random metric patterns obtained by sampling a few points at a time. The order n controls how much structure is visible. At n = 2 we see only a distance histogram; at n = 3 we see random metric triangles; at larger n we see increasingly global finite subspaces. 3.1 Random Distance Matrix Laws (n)
The law ρX introduced above is an intrinsic finite-dimensional invariant of X . It records the distribution of random n-point metric subspaces sampled according to µ, with labels retained only through the ordering of the sampled tuple. Since the samples are i.i.d., (n) ρX is invariant under simultaneous permutations of the n sampled points; we keep the ordered-vector representation because it gives a convenient Euclidean target space for optimal transport. Definition 3.1 (Distance-Matrix Wasserstein) The order-n Distance-Matrix Wasserstein statistic is (n) (n) DMWn,p (X , Y) = Wp (ρX , ρY ), where the Wasserstein distance is computed using ∥ · ∥p,avg on RNn . The order n = 2 recovers a Wasserstein comparison of distance distributions. Higher orders compare random triangles, four-point configurations, and increasingly expressive finite metric subspaces. 3.2 Projection Consistency The distance-matrix laws form a projective family. If 2 ≤ k ≤ n, let Pn→k : RNn → RNk denote the coordinate projection that keeps only the distances among the first k sampled points. Proposition 3.2 (Projection consistency) For every metric measure space X and every 2 ≤ k ≤ n, (n) (k) (Pn→k )# ρX = ρX . Moreover, DMWk,p (X , Y) ≤ DMWn,p (X , Y). The proof of Proposition 3.2 is given in Section C.1. Remark 3.3 (Normalization and hierarchy) For a fixed coordinate projection, the averaged norm can expand by (Nn /Nk )1/p because the denominator changes from Nn to Nk . The monotone hierarchy above avoids this artifact by averaging over all k-submatrices. Thus the averaged normalization remains compatible with a genuine lower-bound hierarchy. 7
Xu and Wu
3.3 Relaxation View Every pointwise coupling π ∈ Π(µ, ν) induces a coupling of distance matrices by sampling (x1 , y1 ), . . . , (xn , yn ) ∼ π independently. Let n o (n) (n) Cn = (TX , TY )# π ⊗n : π ∈ Π(µ, ν) . Then (n)
(n)
Cn ⊆ Π(ρX , ρY ). Thus DMW removes the requirement that a distance-matrix coupling be induced by a single pointwise coupling. Proposition 3.4 (Relaxation characterization) For every n ≥ 2, (n)
(n)
Cn ⊆ Π(ρX , ρY ). Consequently, DMW minimizes the same distance-matrix distortion over a coupling class containing all couplings induced by point-level GW couplings. The proof of Proposition 3.4 is given in Section C.2.
4 Main Theory This section explains why DMW is not merely a descriptive feature of metric data. The argument has two sides. The easy direction is that DMW is a relaxation of GW, hence a lower bound. The harder and more important direction is that this relaxation is not loose in the large-order limit: the gap is controlled by the ordinary empirical Wasserstein error of sampling the underlying spaces. The auxiliary tools needed for this argument are collected in Section A and proved in Section B. The proof of the main approximation theorem proceeds in three steps. First, Theorem 4.1 shows that DMW is a lower bound of GW. Second, the lifting result Lemma B.3 and empirical control lemma Lemma B.4 connect matrix-level couplings back to empirical metric measure spaces. Third, the GW triangle inequality Lemma B.1 and same-space stability estimate Lemma B.2 yield Theorem 4.2. 4.1 Lower Bound Theorem 4.1 (DMW is a GW lower bound) For every n ≥ 2 and p ≥ 1, DMWn,p (X , Y) ≤ GWp (X , Y). The proof of Theorem 4.1 is given in Section C.3. The lower bound is conceptually important but, by itself, would not be enough for a learning method: many cheap summaries are lower bounds. The next result is the point of the construction. It says that if an n-sample empirical measure is already a good Wasserstein approximation of the original space, then an order-n distance-matrix law contains enough information to recover GW up to that sampling error. 8
Distance-Matrix Wasserstein Statistics
4.2 Reverse Approximation Bound Theorem 4.2 (GW approximation by DMW) For every n ≥ 2, GWp (X , Y) ≤ DMWn,p (X , Y) + 2EWp (µ, µ bn ) + 2EWp (ν, νbn ). The proof of Theorem 4.2 is given in Section C.4. Corollary 4.3 (Gap control) 0 ≤ GWp (X , Y) − DMWn,p (X , Y) ≤ 2EWp (µ, µ bn ) + 2EWp (ν, νbn ). The proof of Corollary 4.3 is given in Section C.5. The gap bound gives a clean interpretation of the hierarchy. The order n is not an arbitrary feature dimension; it is the sample size of a random finite subspace. Larger n makes the random subspace a better proxy for the original metric measure space, and the (n) price of this improvement appears later as a harder distribution-estimation problem for ρX . 4.3 Statistical Consequences: Rates of Convergence The preceding theorem reduces the population approximation error of DMW to classical empirical Wasserstein approximation. We record several consequences that will be used to interpret the choice of the order n. Corollary 4.4 (Population consistency) For compact metric measure spaces X and Y, lim DMWn,p (X , Y) = GWp (X , Y).
n→∞
The proof of Corollary 4.4 is given in Section C.6. Corollary 4.5 (Finite-space rate) Assume X and Y are finite with cardinalities mX and mY , and let RX = diam(X) and RY = diam(Y ). Then EWp,X (µ, µ bn ) ≤ RX
√ mX 1/p √ , 2 n
EWp,Y (ν, νbn ) ≤ RY
√ mY 1/p √ . 2 n
Consequently, √ √ mX 1/p mY 1/p √ √ GWp (X , Y) − DMWn,p (X , Y) ≤ 2RX + 2RY . 2 n 2 n In particular, for p = 1 this gives the parametric bound r r mX mY GW1 (X , Y) − DMWn,1 (X , Y) ≤ RX + RY . n n The proof of Corollary 4.5 is given in Section C.7. 9
Xu and Wu
Remark 4.6 (Sharpness of the finite-space bound) For p > 1, the finite-space rate above should be interpreted carefully. The bound is obtained from Wpp ≤ Rp TV and Jensen’s inequality, and it is not a sharp statement about the transportation linear program in every finite metric. However, the exponent n−1/(2p) cannot be uniformly improved for Wp itself without additional structure: on a two-point space with positive masses on both points, Wp (µ, µ bn ) is proportional to |b p − p|1/p , whose expectation is of order n−1/(2p) . The pth p power Wp has the usual parametric behavior. In graph-learning applications we often use p = 1, where the displayed rate for W1 has the familiar n−1/2 form. Remark 4.7 (Typical Euclidean and manifold rates) Under standard compact-support empirical Wasserstein assumptions, the abstract gap bound can be converted into dimensiondependent rates. For example, if µ is supported on a compact subset of Rd of diameter at most R, compact-support bounds for EWpp (µ, µ bn ) imply, after Jensen’s inequality, rates of the form n−1/(2p) , d < 2p, −1/(2p) 1/p EWp (µ, µ bn ) ≤ C n (log n) , d = 2p, −1/d n , d > 2p. for a constant depending on the support and on (d, p); see, for example, Fournier and Guillin (2015). Analogous rates hold on compact manifolds under the usual geometric regularity assumptions. We use these rates as interpretation of the main theorem rather than as new empirical process results.
5 Empirical and Sliced Estimators The preceding section is a population statement: it assumes access to the entire distancematrix laws. In computation we only see sampled distance matrices, and direct Wasserstein estimation in RNn becomes expensive when n grows. This section separates three layers (n) of approximation: empirical estimation of ρX , concentration of the resulting statistic, and sliced projections that replace high-dimensional OT by many one-dimensional OT problems. 5.1 Empirical DMW (n)
(n)
(n)
Let DX,1 , . . . , DX,K be independent samples from ρX and define K
(n)
ρbX,K =
1 X δD(n) . K X,k k=1
The empirical statistic is (n) (n) \ n,p,K = Wp (b DMW ρX,K , ρbY,K ).
Theorem 5.1 (Empirical error decomposition) Let (n)
(n)
(n)
ηX (n, K) = EWp (b ρX,K , ρX ),
(n)
ηY (n, K) = EWp (b ρY,K , ρY ).
Then \ n,p,K − DMWn,p (X , Y) ≤ ηX (n, K) + ηY (n, K). EDMW 10
Distance-Matrix Wasserstein Statistics
The proof of Theorem 5.1 is given in Section C.8. Theorem 5.2 (High-probability empirical DMW bound) Let R = max{diam(X), diam(Y )}. For any δ ∈ (0, 1), with probability at least 1 − δ over the K sampled distance matrices from each space, r 2 1/2−1/p \ n,p,K − DMWn,p (X , Y) ≤ ηX (n, K) + ηY (n, K) + RK DMW log . δ In particular, for p = 1 the stochastic deviation term is r log(2/δ) R . K The proof of Theorem 5.2 is given in Section C.9. Remark 5.3 (Why the bound is sharpest for p = 1) The proof uses a bounded-difference (n) (n) argument for the statistic Wp (b ρX,K , ρbY,K ). Replacing one atom of an empirical measure changes that empirical measure by at most RK −1/p in Wp . This gives the clean parametric concentration rate when p = 1. For p > 1, it is often more natural to concentrate the pth-power objective, as in Theorem 5.4. Theorem 5.4 (Powered concentration) Let R = max{diam(X), diam(Y )}. For any δ ∈ (0, 1), the following bound holds with probability at least 1 − δ: r p p \ n,p,K − EDMW \ n,p,K ≤ Rp log(2/δ) . DMW K The same concentration holds for full empirical DMW at every fixed order n and dimension p \ n,p,K − DMWpn,p , not Nn ; the dependence on dimension enters through the bias term EDMW through the bounded-difference fluctuation. The proof of Theorem 5.4 is given in Section C.10. Remark 5.5 (Why concentrate the powered objective) For p > 1, the natural optimaltransport objective is the transport cost Wpp before taking the pth root. Concentrating this powered objective is standard in computational OT analyses, including entropic OT and Sinkhorn-type losses, because the cost is bounded and changes by only O(K −1 ) when one empirical atom is replaced. The root can always be taken afterward, but doing so generally weakens constants near zero. Theorem 5.6 (Total empirical approximation error) With ηX (n, K) and ηY (n, K) as above, \ n,p,K − GWp (X , Y) ≤ 2EWp,X (µ, µ EDMW bn ) + 2EWp,Y (ν, νbn ) + ηX (n, K) + ηY (n, K). The proof of Theorem 5.6 is given in Section C.11. 11
Xu and Wu
Corollary 5.7 (High-probability total approximation) Under the assumptions of Theorem 5.2, with probability at least 1 − δ, \ n,p,K − GWp (X , Y) ≤ 2EWp,X (µ, µ DMW bn ) + 2EWp,Y (ν, νbn ) r + ηX (n, K) + ηY (n, K) + RK
1/2−1/p
2 log . δ
The proof of Corollary 5.7 is given in Section C.12. ′ Proposition 5.8 (Permutation validity) Let Z1 , . . . , Zm and Z1′ , . . . , Zm ′ be two samples of metric measure spaces. Let T be any statistic computed from the two labeled groups, such as DMW, sliced DMW, multi-scale DMW, or a DMW kernel MMD. Under the null hypothesis that the two groups are exchangeable, the exact permutation test that rejects for large values of T has type-I error at most α.
The proof of Proposition 5.8 is given in Section C.13. Corollary 5.9 (Direct DMW estimation in ambient dimension) Let RX = diam(X) (n) (n) and RY = diam(Y ), and set R = max{RX , RY }. Since ρX and ρY are supported in [0, RX ]Nn and [0, RY ]Nn , standard empirical Wasserstein bounds in RNn imply, for a constant C = C(p, Nn , R), −1/(2p) , Nn < 2p, K −1/(2p) 1/p ηX (n, K) + ηY (n, K) ≤ C K (log K) , Nn = 2p, −1/Nn K , Nn > 2p. Thus direct empirical DMW estimation suffers from the ambient dimension Nn = n(n−1)/2. The proof of Corollary 5.9 is given in Section C.14. Theorem 5.10 (Intrinsic-dimensional estimation via covering numbers) Assume that the supports of µ and ν satisfy polynomial covering bounds N (supp µ, dX , ε) ≤ AX ε−sX ,
N (supp ν, dY , ε) ≤ AY ε−sY
for 0 < ε ≤ 1. Then, for every fixed order n, the supports of the matrix laws satisfy nsX 2 (n) n N (supp ρX , ∥ · ∥p,avg , ε) ≤ AX , ε and analogously for Y . Consequently, there are constants CX (n, p) and CY (n, p) such that ηX (n, K) ≤ CX (n, p) ψnsX ,p (K), where
ηY (n, K) ≤ CY (n, p) ψnsY ,p (K),
−1/(2p) , K ψs,p (K) = K −1/(2p) (log K)1/p , −1/s K ,
s < 2p, s = 2p, s > 2p.
Consequently, the matrix-law estimation term in Theorem 5.6 scales with the intrinsic product dimensions nsX and nsY , rather than the ambient matrix dimension Nn = n(n − 1)/2. 12
Distance-Matrix Wasserstein Statistics
The proof of Theorem 5.10 is given in Section C.22. Remark 5.11 (Why intrinsic dimension matters) Theorem 5.10 is the statistical reason DMW can remain meaningful when the distance-matrix vector has many coordinates. The ambient space of n-point distance matrices has dimension Nn = n2 , but those matrices are generated by n points moving in the original space. If the original data live on a d-dimensional geometric object, the image of the distance-matrix map has effective dimension at most order nd, not n2 . This does not remove the approximation–estimation tradeoff, but it explains why empirical DMW can be far less pessimistic than a naive high-dimensional Wasserstein bound suggests. Corollary 5.12 (Compact manifold sample complexity) If X and Y are compact Riemannian manifolds of dimensions sX and sY with bounded geometry, for example with uniformly bounded diameter, Ricci curvature bounded below, and injectivity radius bounded away from zero, then the covering assumption of Theorem 5.10 holds. Hence empirical DMW has the intrinsic-dimensional rate displayed in Theorem 5.10. The proof of Corollary 5.12 is given in Section C.23. Remark 5.13 (Scope of the intrinsic-dimensional rate) Theorem 5.10 should be read as a fixed-order estimation result. It mitigates the ambient dimension Nn for structured data, but the effective dimension still grows linearly in n. This is exactly the approximation– estimation tradeoff: increasing n improves the population approximation to GW, while es(n) timating ρX becomes harder. 5.2 Sliced DMW Sliced DMW addresses the computational bottleneck directly. Rather than transport two empirical distributions in the full distance-matrix space, we project each distance-matrix vector onto many random one-dimensional directions and compute one-dimensional Wasserstein distances by sorting. The following normalization only ensures that every projection is 1-Lipschitz for the averaged matrix norm. Let q be the Holder conjugate of p, with the usual convention q = ∞ when p = 1. The dual norm of ∥ · ∥p,avg is ∥θ∥q,dual = Nn1/p ∥θ∥q . Let Θn,p = {θ ∈ RNn : ∥θ∥q,dual ≤ 1} and let σ be a probability distribution on Θn,p . For θ ∈ Θn,p , define Pθ (a) = ⟨θ, a⟩. Then Pθ is 1-Lipschitz from (RNn , ∥ · ∥p,avg ) to R. Define h i (n) (n) 1/p SDMWn,p (X , Y) = Eθ∼σ Wpp (Pθ )# ρX , (Pθ )# ρY . The empirical Monte Carlo version with L directions θ1 , . . . , θL ∼ σ is " #1/p L X 1 (n) (n) \ n,p,K,L = SDMW Wpp (Pθℓ )# ρbX,K , (Pθℓ )# ρbY,K L ℓ=1
13
Xu and Wu
. Each one-dimensional Wasserstein distance is computed by sorting. Proposition 5.14 (Sliced lower hierarchy) For every n ≥ 2, SDMWn,p (X , Y) ≤ DMWn,p (X , Y) ≤ GWp (X , Y). The proof of Proposition 5.14 is given in Section C.15. Proposition 5.15 (Cramer–Wold expressivity of sliced distance laws) Fix n ≥ 2. If (n) (n) (Pθ )# ρX = (Pθ )# ρY for every θ ∈ RNn , (n)
(n)
then ρX = ρY . Consequently, if this equality of all one-dimensional projected laws holds for every n ≥ 2, then X and Y are measure-preserving isometric. The proof of Proposition 5.15 is given in Section C.24. Corollary 5.16 (Sliced hierarchy is complete with all directions) Assume that the slicing law σ has full support on Θn,p for every n. If SDMWn,p (X , Y) = 0
for every n ≥ 2,
then X and Y are measure-preserving isometric. The proof of Corollary 5.16 is given in Section C.25. Remark 5.17 (Information loss from finite slicing) Proposition 5.15 is an all-directions (n) statement. A finite Monte Carlo set of directions does not determine ρX and should be interpreted as a scalable estimator of the sliced population statistic, not as a complete invariant. If the direction law has full support, equality of projected laws for σ-almost every direction also implies equality of the matrix laws under standard continuity of characteristic functions, but finite L retains only finitely many one-dimensional summaries. Theorem 5.18 (Finite-direction sliced error) Let R = max{diam(X), diam(Y )} and define the empirical sliced population objective (n) (n) SDMWpn,p,K = Eθ∼σ Wpp (Pθ )# ρbX,K , (Pθ )# ρbY,K .
Conditional on the empirical distance matrices, p Rp \ n,p,K,L − SDMWp √ . ≤ Eθ1 ,...,θL SDMW n,p,K 2 L
Consequently the same bound holds after taking expectation over the sampled distance matrices. The proof of Theorem 5.18 is given in Section C.16. 14
Distance-Matrix Wasserstein Statistics
Corollary 5.19 (Finite-direction separation) With the notation of Theorem 5.18, suppose the population sliced signal satisfies, for some ε > 0, SDMWpn,p (X , Y) ≥ εp . For γ, δ ∈ (0, 1), if L≥
R2p 1 log , 2γ 2 ε2p δ
then, with probability at least 1 − δ over the slicing directions, L
1X p (n) (n) Wp (Pθℓ )# ρX , (Pθℓ )# ρY ≥ (1 − γ)εp . L ℓ=1
The proof of Corollary 5.19 is given in Section C.17. Proposition 5.20 (Sliced empirical distribution error) With R = max{diam(X), diam(Y )}, ESDMWpn,p,K − SDMWpn,p (X , Y) ≤ pRp−1 ηX (n, K) + ηY (n, K) . Moreover, for each fixed θ, the one-dimensional projected empirical laws satisfy standard one-dimensional rates; in particular, for p = 1 and compact support of length at most R, (n)
(n)
EW1 (Pθ )# ρbX,K , (Pθ )# ρX
≤ CRK −1/2
for a universal constant C. The proof of Proposition 5.20 is given in Section C.18. Theorem 5.21 (Sliced W1 estimation) Assume p = 1, and set R = max{diam(X), diam(Y )}. Suppose that for every θ ∈ Θn,1 , (n)
(n)
EW1 (Pθ )# ρbX,K , (Pθ )# ρX
≤ C1 RK −1/2 ,
with the analogous bound for Y . Then R \ n,1,K,L (X , Y) − SDMWn,1 (X , Y) ≤ 2C1 RK −1/2 + √ E SDMW . 2 L The proof of Theorem 5.21 is given in Section C.19. 5.3 Multi-Scale DMW No single order n is uniformly best. Small orders are stable and inexpensive but see only local distance statistics; large orders capture more global structure but are harder to estimate. Multi-scale DMW turns this into a controlled design choice by averaging several orders. The theorems below make explicit what this averaging can and cannot do: it preserves lower-bound and approximation guarantees, but consistency still requires the weight distribution to move toward larger orders in the population limit. 15
Xu and Wu
P For a finite index set N ⊂ {2, 3, . . .} and weights αn ≥ 0 with n∈N αn = 1, define the multi-scale sliced statistic X MSDMWp (X , Y) = αn SDMWn,p (X , Y). n∈N
We also write MDMWp (X , Y) =
X
αn DMWn,p (X , Y)
n∈N
for the corresponding unsliced multi-scale DMW. Small orders capture local distance statistics, while larger orders capture more global metric skeletons. Theorem 5.22 (Multi-scale lower bound and approximation) Let N be finite and P let αn ≥ 0 satisfy n∈N αn = 1. Then MSDMWp (X , Y) ≤ MDMWp (X , Y) ≤ GWp (X , Y). Moreover, 0 ≤ GWp (X , Y) − MDMWp (X , Y) ≤
X
αn [2EWp,X (µ, µ bn ) + 2EWp,Y (ν, νbn )] .
n∈N
Consequently, for any sequence of weights αm,n whose mass escapes to high orders, in the sense that X αm,n → 0 for every fixed N, n<N
the corresponding unsliced multi-scale statistics MDMWp,m = GWp on compact metric measure spaces.
P
n αm,n DMWn,p converge to
The proof of Theorem 5.22 is given in Section C.28. Remark 5.23 (Choosing multi-scale weights) Theorem 5.22 shows that multi-scale DMW inherits the single-order GW approximation guarantee by averaging the single-order errors. It does not imply that an arbitrary fixed finite set of small orders is consistent for GW. In practice, weights such as αn ∝ n−1 over a bounded range trade local sensitivity against global structure; for population consistency, the range must grow and the weights must eventually place most mass on large orders. Thus multi-scale DMW is best viewed as a bias–variance design choice rather than a universally faster approximation theorem. Proposition 5.24 (Oracle multi-scale weights) Fix a finite scale set N and suppose that, for each n ∈ N , an upper bound Bn is available for the total order-n error, including population approximation and empirical estimation terms. Then the linear upper bound X X αn Bn , αn ≥ 0, αn = 1, n∈N
n∈N
is minimized by placing all mass on any scale n⋆ ∈ arg minn∈N Bn .
If one imposes an entropy regularization parameter τ > 0 to avoid a degenerate single-scale choice, the optimizer of X X αn Bn + τ αn log αn n∈N
n∈N
16
Distance-Matrix Wasserstein Statistics
is
exp(−Bn /τ ) . m∈N exp(−Bm /τ )
αn⋆ = P
The proof of Proposition 5.24 is given in Section C.29. Theorem 5.25 (Sliced multi-scale closed loop) Assume that for each n the slicing law has full support on Θn,p . Let Xm and X be compact metric measure spaces with uniformly bounded diameters. If SDMWn,p (Xm , X ) → 0
for every fixed n ≥ 2,
then Xm → X in the Gromov-weak topology. If, in addition, the family is precompact in GWp and the limiting Gromov-weak topology agrees with GWp on this family, then GWp (Xm , X ) → 0. In particular, an all-directions, all-orders sliced hierarchy is topologically faithful, although a fixed finite collection of directions and orders is not. The proof of Theorem 5.25 is given in Section C.26. Remark 5.26 (No unconditional numeric equality with GW) Theorem 5.25 is a consistency statement, not a claim that a finite sliced lower bound equals GW. Since SDMWn,p ≤ DMWn,p , sliced multi-scale DMW is generally a further relaxation. Quantitative inequalities of the form DMWn,p ≤ Φn (SDMWn,p ) require inverse Radon or Cramer–Wold stability constants in dimension Nn , which are typically poor. We therefore use sliced DMW as a scalable faithful statistic in the limit of all directions and all orders. We do not treat it as an unconditional finite-sample numerical surrogate equal to GW. 5.4 Kernelization For graph classification and two-sample testing, one may transform a DMW-type dissimilarity into a kernel. A metric alone is not enough to guarantee that exp(−λd) is positive definite. The correct sufficient condition is conditional negative definiteness. Proposition 5.27 (Positive definite sliced DMW kernels) For p = 1, the order-n sliced statistic SDMWn,1 is conditionally negative definite on compact metric measure spaces. Consequently, for every λ > 0, kn (X , Y) = exp{−λSDMWn,1 (X , Y)} is a positive definite kernel. Moreover, any nonnegative weighted multi-scale sum X MSDMW1 (X , Y) = αn SDMWn,1 (X , Y) n∈N
is conditionally negative definite, and therefore kN (X , Y) = exp{−λMSDMW1 (X , Y)} is positive definite. The positive-definiteness conclusion holds for any fixed slicing law with finite first moment on directions; support in Θn,1 is needed for the lower-bound comparison with DMW, not for conditional negative definiteness. 17
Xu and Wu
Method
Leading cost after distance queries are available
Exact empirical DMW Interior-point empirical DMW Sinkhorn empirical DMW Sliced DMW Multi-scale sliced DMW
O(Kn2 + K 2 n2 + OTK ) e 3 )) O(Kn2 + K 2 n2 + O(K 2 2 2 O(Kn + K n + T K 2 ) 2 + LK log K) O(Kn P 2 n∈N O(Kn + LK log K)
Table 1: Representative computational costs. Here K is the number of sampled distance matrices, n is the tuple size, L is the number of slicing directions, T is the number of Sinkhorn iterations, and OTK denotes the cost of the chosen exact balanced OT solver on two K-point empirical measures after the ground-cost matrix has been formed.
The proof of Proposition 5.27 is given in Section C.27. Remark 5.28 (General DMW kernels may be indefinite) Proposition 5.27 is a sufficient condition, not a theorem for every DMW variant. For general p, for unsliced highdimensional DMW, or for nonlinear transformations of DMW, we do not claim positive definiteness unless conditional negative definiteness is verified. In experiments, such similarities should be treated either as dissimilarities, as precomputed indefinite kernels with the usual spectral corrections, or calibrated by permutation procedures rather than by RKHSMMD theory.
6 Algorithms The theoretical sections define a hierarchy of population statistics and prove that they approximate GW. This section explains how the statistics are computed in the regime where GW itself becomes expensive. The key algorithmic idea is a decoupling: GW optimizes over a coupling on the original point sets, so its cost is tied to the number of graph nodes or point-cloud samples. DMW first samples K random n-point subspaces and transports their distance-matrix laws. The scalable variant then slices these matrix laws, replacing a highdimensional K-atom OT problem by L one-dimensional sorting problems. Thus the main computational knobs are the tuple order n, the number of sampled matrices K, and the number of directions L, rather than the full size of the original object once distances can be queried. 6.1 Computational Complexity and Bottlenecks For one sampled n-tuple, forming the full distance vector costs O(n2 ) once pairwise distances can be queried in O(1) time. Thus generating K distance matrices costs O(Kn2 ) per space, or the same order after absorbing the factor of two. If the original graph metric is not precomputed, one must also account for the cost of shortest-path preprocessing; this cost is shared by all DMW variants and is separate from the matrix-law comparison itself. 18
Distance-Matrix Wasserstein Statistics
For full empirical DMW, the two empirical matrix laws have K atoms each in RNn . Building the K × K ground-cost matrix costs O(K 2 Nn ) if distances between matrix vectors are computed naively; since Nn = O(n2 ), this is O(K 2 n2 ). The first practical bottleneck is therefore memory: the full ground-cost matrix already requires O(K 2 ) storage before an OT solver is called. The remaining cost depends on the exact OT solver and is denoted by OTK in Table 1. Following standard complexity bounds for interior-point methods for linear programming, a balanced transport problem on two K-point empirical measures e 3 ) up to accuracy and bitadmits polynomial worst-case arithmetic bounds of order O(K complexity factors. Network simplex implementations are often fast in practice but do not provide the same clean worst-case polynomial guarantee, so we do not use network-simplex complexity as a theorem-level statement. Entropic empirical DMW replaces exact OT by Sinkhorn iterations. Once the cost matrix has been formed, T Sinkhorn iterations cost O(T K 2 ) arithmetic operations and O(K 2 ) memory. In this paper, Sinkhorn is used only as a numerical approximation to empirical DMW. Standard entropic OT theory gives convergence of Sinkhorn iterates to the unique entropically regularized coupling for fixed regularization parameter (Cuturi, 2013; Peyré and Cuturi, 2019), and the regularized objective converges to the unregularized empirical OT value as the regularization parameter tends to zero under the usual compact-cost assumptions. We do not analyze gradient-flow training dynamics here, because DMW is used as a statistic or kernel rather than as an end-to-end learned loss in the main experiments. 6.2 Efficient Implementation via Slicing and Streaming Sliced DMW avoids the K × K high-dimensional OT solve. After the O(Kn2 ) sampling step, each slicing direction requires projections of the K distance vectors and two sorts. If projected values are formed from stored distance vectors, this costs O(KNn + K log K) per direction. If projections are accumulated while generating distance entries, the arithmetic can be organized as O(Kn2 + LK log K) up to the cost of drawing and applying directions. This realizes the statistical estimator of Theorems 5.18 and 5.21 through one-dimensional Wasserstein computations. The implementation is naturally vectorized. With distance vectors arranged as a matrix A ∈ RK×Nn and directions as Θ ∈ RL×Nn , all projections are computed as the matrix product AΘ⊤ . The L sorting problems are independent and can be performed in parallel across directions. Modern array frameworks such as PyTorch, JAX, or GPU-backed NumPy libraries can therefore batch the projection step and parallelize sorting across L. This is not needed for the moderate-scale experiments in this paper, but it is important for the intended large-scale setting: the sliced estimator exposes parallelism that the K × K OT solve of full DMW does not. The memory footprint also admits a streaming implementation. Storing all distance vectors costs O(KNn ), which is usually acceptable for the tuple sizes used in the experiments. If memory is tight, one can process one direction, or a small batch of directions, at a time: sample an n-tuple, compute its distance entries, accumulate its projected values for the active directions, and discard the tuple. This reduces the working memory to O(K) per active direction, plus the memory needed for distance queries. Streaming does not change the estimator; it changes only how the same projected samples are materialized. 19
Xu and Wu
Algorithm 1 Sliced Distance-Matrix Wasserstein Require: Metric measure spaces X = (X, dX , µ) and Y = (Y, dY , ν); distance oracles or precomputed distance matrices for dX , dY ; tuple size n; number of tuple samples K; number of slicing directions L; exponent p; direction law σ on Θn,p . \ n,p,K,L (X , Y). Ensure: Empirical sliced statistic SDMW 1: // Step 1: Sample finite metric subspaces 2: for k = 1 to K do i.i.d. 3: Draw xk,1 , . . . , xk,n ∼ µ and set Ak ← dX (xk,i , xk,j ) 1≤i<j≤n ∈ RNn . 4:
i.i.d.
Draw yk,1 , . . . , yk,n ∼ ν and set Bk ← dY (yk,i , yk,j ) 1≤i<j≤n ∈ RNn .
5: end for 6: // Step 2: Project distance matrices to one dimension 7: for ℓ = 1 to L do 8: 9: 10: 11: 12:
Draw a slicing direction θℓ ∼ σ from the dual unit set Θn,p . for k = 1 to K do aℓk ← ⟨θℓ , Ak ⟩, bℓk ← ⟨θℓ , Bk ⟩. end for Sort the projected samples: aℓ(1) ≤ · · · ≤ aℓ(K) ,
13:
bℓ(1) ≤ · · · ≤ bℓ(K) .
Compute the one-dimensional transport cost K
1 X cℓ ← |aℓ(k) − bℓ(k) |p . K k=1
14: end for 15: // Step 3: Average over slicing directions
\ n,p,K,L ← 16: SDMW
P L 1 L
ℓ=1 cℓ
1/p
.
\ n,p,K,L . 17: return SDMW
6.3 Practical Heuristics: Sampling and Budgets The parameter K should be interpreted as a Monte Carlo accuracy budget, not as the number of nodes in the original graph or point cloud. The population theory fixes the underlying metric measure space and lets K control empirical estimation of the distancematrix law. In the scalability experiments we deliberately keep K fixed while increasing the graph size in order to isolate the cost of avoiding node-level GW couplings. If one requires the empirical DMW estimation error to decrease with the original dataset size, then K must 20
Distance-Matrix Wasserstein Statistics
also be increased according to the desired statistical accuracy; in that regime the LK log K sorting term can become the new bottleneck. This is the intended tradeoff: DMW replaces an all-node coupling problem by an explicit, user-controlled sampling budget. For graphs, there is a second budget: how pairwise distances are made available. On small and medium graphs we precompute all-pairs shortest-path distances. On very large graphs, one may instead sample a node subset, use landmark distances, or compute shortest paths on demand. Such procedures should be viewed as approximations to the original metric measure space before DMW is applied. They introduce a separate subgraph or distance-oracle error that is not part of the population DMW theorem. This distinction is important in the experiments: node-budget truncation is an implementation device for large social graphs, while the formal DMW estimator assumes sampling from the original measure and exact access to the chosen metric. The tuple order n and scale set N are statistical choices. Larger n improves the population approximation to GW, but Fig. 3 shows that finite-K variance can dominate at high order. In practice, we therefore use modest orders such as {2, 3, 4, 6} and tune K, L, and the scale weights inside training folds when labels are available. This is also why the multi-scale version is used in the large experiments: it turns the bias–variance behavior predicted by the theory into a controlled algorithmic design.
7 Expressivity The approximation theorem says that increasing n recovers GW at the population level. It is equally important to understand what is lost at any fixed order. Fixed-order DMW is a compressed statistic: it records the law of random n-point metric patterns, not the whole space. This limitation is already visible for n = 2. Proposition 7.1 (A fixed-order counterexample) There exist two finite metric measure spaces X and Y with uniform measures such that (2)
(2)
ρX = ρY . Nevertheless, X and Y are not measure-preserving isometric. The example used in the proof has four points in each space. In one space, the pairs at distance two form a matching, equivalently the opposite pairs of a four-cycle. In the other, the pairs at distance two share a common endpoint. Both spaces therefore have the same unordered distance histogram, but their distance-two graphs have different degree patterns. This makes the example a minimal illustration of why order-2 DMW cannot be complete. The proof of Proposition 7.1 is given in Section C.20. Proposition 7.2 (Completeness of the full hierarchy) Let X and Y be compact metric measure spaces. If (n) (n) ρX = ρY for every n ≥ 2, then X and Y are measure-preserving isometric. The proof of Proposition 7.2 is given in Section C.21. Thus fixed-order DMW trades completeness for statistical and computational tractability, while the full distance-matrix hierarchy is complete. 21
Xu and Wu
Y
X
distance-two graph distance-two graph is a matching shares a vertex Both spaces have the same order-2 distance law; the distance-two degree patterns differ.
Figure 2: Distance-two graphs for the fixed-order counterexample. Both spaces have exactly two unordered pairs at distance two and four unordered pairs at distance one, hence the same order-2 distance distribution. The degree patterns of the distancetwo graphs differ, so the spaces are not measure-preserving isometric.
8 Experiments The experiments are designed as empirical tests of the theory rather than as a single benchmark table. We evaluate four claims: increasing the order n improves structural expressivity but increases estimation difficulty; random slicing gives a controlled accuracy–runtime knob; multi-scale sliced DMW is computationally scalable; and the resulting positivedefinite p = 1 kernel is useful, but not universally dominant, for graph learning and twosample testing. The reproducibility repository is available at https://github.com/Ao-Xu/ dmw-reproducibility. 8.1 Protocol and Reproducibility Graphs are represented as metric measure spaces using normalized shortest-path distance and the uniform node measure. For very large social graphs, the implementation caps the metric computation at a fixed node budget when constructing DMW features; this is a finite-subgraph approximation used by the scalable algorithm, and the full graph labels are still used in the benchmark. This engineering step is separate from the population theory: theorems such as Theorem 4.2 assume i.i.d. sampling from the original measure, whereas node-budget truncation first replaces the original graph by a sampled finite subgraph and therefore introduces an additional approximation error. We use the truncation only as a large-graph implementation device, not as part of the formal DMW definition. Sliced DMW samples tuples with replacement, uses Euclidean-normalized Gaussian directions, and computes one-dimensional Wasserstein distances by sorting. The dual-unit normalization in Proposition 5.14 is used for theorem-level lower bounds; Euclidean normalization changes only the scale of the sliced statistic and is absorbed by kernel bandwidth selection. Unless otherwise stated, we use p = 1 so that the exponential sliced multi-scale kernel is positive definite by Proposition 5.27. For graph classification, all SVM hyperparameters are selected inside the training folds. Small and medium TU datasets use 10-fold nested cross-validation. For the largest datasets, 22
Distance-Matrix Wasserstein Statistics
Empirical validation of the approximation--estimation tradeoff
(a) Population proxy improves with order
(b) Estimation variability grows 0.035
0.028
Std. over repetitions
Approximation gap
0.030
0.026 0.024 0.022 0.020 0.018
0.030 0.025 0.020 0.015 0.010
0.016
(c) Finite-sample tradeoff
(d) Estimates with uncertainty 0.14 0.13
SDMW estimate
Total empirical error
0.035
0.030
0.025
0.020
0.12 0.11 0.10 0.09 0.08 0.07
2
4
6
8
10
12
Distance-matrix order n
2
4
6
8
10
12
Distance-matrix order n
Figure 3: Empirical validation of the approximation–estimation tradeoff. Higher order improves the population proxy but increases finite-sample variability; the total empirical error is minimized at an intermediate order.
NCI1 and REDDIT-BINARY, we use the full dataset with a reduced inner validation grid and 5 outer folds; this keeps the benchmark computationally feasible while preserving an honest train-test split. The baselines are shortest-path histograms, degree histograms, Weisfeiler–Lehman subtree features, graphlet/motif features, NetLSD spectral features, and entropic GW kernels on small graph datasets where pairwise GW computation is feasible. 8.2 Theory Validation: The Approximation–Estimation Tradeoff Theorem 4.2 and Corollary 4.3 predict that the population DMW approximation improves with order, while Corollary 5.9 and Theorem 5.10 predict that the empirical matrix law becomes harder to estimate at larger order. We test this on noisy circle-versus-ellipse metric spaces. A high-order, high-sample sliced statistic is used as a reference, and lower-order statistics are repeatedly estimated with a fixed Monte Carlo budget. Figure 3 is the main sanity check for the paper’s statistical story. It shows why a single large order is not automatically best in finite samples. The multi-scale statistic is 23
Xu and Wu
Finite-direction concentration of sliced DMW (a) Separation probability
(b) Monte Carlo variability
1.0
0.8
Sliced statistic
Within 20% of high-L reference
0.07
0.6
0.4
0.06
0.05
0.04
0.2 Δ = 0.03 Δ = 0.05
0.0 2
1
2
2
Δ = 0.08 Δ = 0.12
2
3
2
4
2
0.03 5
2
6
2
7
2
Number of directions L
1
2
2
2
3
2
4
2
5
2
6
2
7
Number of directions L
Figure 4: Finite-direction concentration. Increasing L improves the probability of approximating the high-direction sliced statistic and reduces Monte Carlo variability.
designed precisely for this regime: it aggregates several orders rather than betting the entire estimator on one high-dimensional distance-matrix law. Quantitatively, the reference gap decreases from about 3.1 · 10−2 at order n = 2 to about 1.7 · 10−2 at order n = 12 in this run. This is consistent with the population approximation theorem. The empirical standard deviation, however, remains visible at larger orders under a fixed sampling budget, which is the finite-sample effect predicted by the matrix-law estimation bounds. Thus the figure does not merely illustrate a heuristic; it checks the two competing terms in the theory. 8.3 Finite-Direction Slicing Corollary 5.19 predicts that the number of slicing directions L controls the probability of retaining a population separation signal. We compare circle-versus-ellipse metric spaces with different signal strengths and measure how often an L-direction estimate lies within 20% of a high-L reference. The monotone stabilization in Fig. 4 is consistent with the L−1/2 concentration term. In practice, L is therefore an explicit accuracy knob: small values can be enough for strong signals, while weaker signals require more directions. 8.4 Scalability and Parameter Bottlenecks We compare exact POT GW, entropic GW, full empirical DMW, and sliced DMW on stochastic block model graph metrics. Exact GW is reported only where the solver remains reliable; missing entries are not extrapolated. The DMW methods use fixed Monte Carlo budgets, so their runtime is controlled by tuple sampling and sorting rather than by optimizing a node-level coupling. The missing runtime entries have a specific interpretation. They do not indicate failed DMW runs. They mark GW baselines that were intentionally 24
Distance-Matrix Wasserstein Statistics
Scalability: DMW avoids node-level GW optimization Full-DMW POT-GW POT-entropic-GW Sliced-DMW
Wall-clock time (s)
10
10
−1
−2
10
2
10
3
Number of graph nodes
Figure 5: Runtime scaling on SBM graph metrics. Sliced DMW avoids the node-level GW optimization bottleneck and remains stable up to 1000 nodes with fixed Monte Carlo budgets.
Runtime vs n
10
−2
10
Seconds
10
Runtime vs L
Seconds
Seconds
Runtime vs K
−2
−2
9 × 10 8 × 10
10
−3
6 × 10 2
5
2
7
2
9
7 × 10
−3
2
4
2
K
6
L
2
8
−3
−3
−3
2 × 10
0
0
3 × 104 × 10
0
6 × 10
0
10
1
Order n
Figure 6: Runtime as a function of the core Monte Carlo parameters. The empirical behavior matches the implementation cost O(Kn2 + LK log K) once graph distances are available.
not run beyond the regime where a node-level GW solver is a meaningful scalable comparator. Exact GW is omitted beyond 120 nodes, and entropic GW is omitted at 1000 nodes, because those baselines require optimizing over couplings on the full node sets. We leave these cells blank rather than extrapolate runtimes from smaller graphs. 25
Xu and Wu
Dataset
MS-SDMW
Entropic-GW
Shortest-path
Degree
WL subtree
Graphlet
NetLSD
BZR COX2 ENZYMES IMDB-BINARY IMDB-MULTI MUTAG NCI1 PROTEINS PTC MR REDDIT-BINARY
0.785 ± 0.019 0.782 ± 0.009 0.205 ± 0.082 0.583 ± 0.048 0.391 ± 0.022 0.767 ± 0.092 0.538 ± 0.022 0.708 ± 0.046 0.567 ± 0.039 0.779 ± 0.013
0.788 ± 0.011 0.782 ± 0.009 – – – 0.686 ± 0.076 – – 0.523 ± 0.037 –
0.802 ± 0.042 0.782 ± 0.009 0.270 ± 0.062 0.593 ± 0.054 0.399 ± 0.024 0.777 ± 0.076 0.605 ± 0.013 0.715 ± 0.025 0.578 ± 0.039 0.772 ± 0.017
0.790 ± 0.025 0.782 ± 0.009 0.265 ± 0.061 0.578 ± 0.042 0.415 ± 0.022 0.803 ± 0.092 0.598 ± 0.010 0.668 ± 0.029 0.570 ± 0.032 0.744 ± 0.023
0.887 ± 0.055 0.844 ± 0.071 0.560 ± 0.080 0.729 ± 0.041 0.510 ± 0.043 0.846 ± 0.059 0.854 ± 0.011 0.758 ± 0.033 0.663 ± 0.071 0.825 ± 0.022
0.872 ± 0.058 0.790 ± 0.020 0.398 ± 0.061 0.700 ± 0.039 0.465 ± 0.025 0.884 ± 0.081 0.654 ± 0.019 0.749 ± 0.044 0.576 ± 0.064 0.886 ± 0.015
0.849 ± 0.038 0.773 ± 0.023 0.355 ± 0.067 0.702 ± 0.032 0.470 ± 0.047 0.857 ± 0.056 0.676 ± 0.016 0.739 ± 0.038 0.590 ± 0.049 0.853 ± 0.010
Table 2: TU graph classification. Entries report mean accuracy ± standard deviation. MSSDMW is competitive with metric baselines and entropic GW on several datasets, but WL and graphlet features dominate on label-rich combinatorial benchmarks.
Figures 5 and 6 separate two computational effects. First, replacing node-level GW by distance-matrix sampling removes the all-node coupling bottleneck. Second, the remaining costs are transparent: K controls how many distance matrices are sampled, L controls how many one-dimensional projections are sorted, and n controls the cost of forming each tuple distance vector. 8.5 TU Graph Classification We evaluate on ten TU datasets covering molecular, protein, and social graphs. The suite consists of MUTAG, PTC MR, BZR, COX2, PROTEINS, ENZYMES, IMDB-BINARY, IMDB-MULTI, NCI1, and REDDIT-BINARY. The datasets range from 188 to 4110 graphs, with average graph sizes from about 13 to 430 nodes. The MS-SDMW kernel uses orders {2, 3, 4, 6} or {2, 3, 4, 6, 8} depending on the dataset size, with the budgets recorded in the experiment script. The dashes in the Entropic-GW column have the same computational meaning as the missing GW entries in the runtime experiment. A GW kernel requires computing a GW distance for every pair of graphs in the dataset. For NCI1, for example, this would require more than eight million pairwise GW computations; for REDDIT-BINARY, individual graphs can contain thousands of nodes. We therefore report Entropic-GW only on the small molecular datasets where a pairwise GW kernel can be computed without changing the benchmark protocol. The omission is part of the scalability comparison, not an unavailable label or a failed run. Table 2 and Fig. 7 give a deliberately conservative result. MS-SDMW is not state-of-theart on motif-driven graph classification. Its strongest empirical role is more specific: it is a scalable, GW-grounded metric kernel. This distinction matters for the paper’s claim. The method is designed for intrinsic metric comparison, not for replacing WL on datasets where discrete labels and local refinement features explain most of the task. This distinction also explains the pattern across datasets. On PROTEINS and REDDIT-BINARY, MSSDMW is competitive with shortest-path and degree-based metric baselines. On ENZYMES and NCI1, the gap to WL is large, which indicates that node labels and local discrete refinement carry information that is not present in the metric-only representation used here. These results support the theory where it applies, namely intrinsic metric comparison, 26
Distance-Matrix Wasserstein Statistics
MS-SDMW Entropic-GW
NetLSD
0.83 0.78
0.77 0.74
0.65 0.68
0.70 0.70
0.73
0.76 0.75 0.74
0.72
0.71
0.58
0.59 0.58
0.54
0.56
0.61 0.60
0.67
0.66 0.58 0.59
0.58 0.57
0.51
0.52
0.57
0.6
0.85
0.85
0.84 0.79 0.77
0.78 0.78 0.78 0.78
0.79 0.79 0.80 0.79
0.89
0.89 0.87 0.85
0.86
0.88 0.77 0.69
Nested-CV accuracy
0.78 0.80
0.85
0.9
0.7
WL subtree Graphlet
TU graph classification: metric DMW kernel versus structural baselines
1.0
0.8
Shortest-path Degree
0.40 0.41
0.35
0.4
0.39
0.40
0.47 0.47
0.5
G TA MU
MR C_
PT
R
BZ
CO
X2
NS
ME
EI
PR
OT
EN
S
NA
ZY
I -B DB
RY B MD
U -M
I LT
NC
I1
I
IM
RE
NA BI TDI
RY
D
Figure 7: Expanded TU benchmark across ten datasets. The figure visualizes the full table and makes the main limitation clear: DMW captures metric structure, whereas WL and graphlet kernels directly exploit discrete motifs and node labels.
Accuracy--runtime profile on TU benchmarks 0.9
Nested-CV accuracy
0.8 0.7 0.6 0.5 method MS-SDMW WL subtree Graphlet NetLSD dataset MUTAG PTC_MR
0.4 0.3 0.2 10
−2
10
−1
10
0
10
1
BZR COX2 PROTEINS ENZYMES IMDB-BINARY IMDB-MULTI NCI1 REDDIT-BINARY
10
2
Kernel/feature construction time (s, log scale)
Figure 8: Accuracy–runtime profile on TU datasets. MS-SDMW sits in the unsupervised metric-kernel regime: cheaper than pairwise GW kernels and often competitive with metric baselines, but not a universal replacement for specialized graph kernels.
while identifying fused or label-aware DMW as the right next extension for label-rich graph learning. 27
Xu and Wu
Two-sample power of the MS-SDMW kernel 1.0
Rejection probability
0.8
0.6
Δ = 0.00 Δ = 0.04 Δ = 0.08
Δ = 0.12 Δ = 0.16
Δ = 0.24 nominal 0.05
0.4
0.2
0.0 20
40
60
80
100
Metric spaces per group
Figure 9: Two-sample testing power for the MS-SDMW kernel. The type-I error remains near the nominal 0.05 level under the null, and power increases with both sample size and geometric signal.
8.6 Two-Sample Testing We use the positive-definite MS-SDMW kernel in an MMD two-sample test. Samples are noisy point clouds drawn from a circle or an ellipse with eccentricity shift ∆. The null case is ∆ = 0. P-values are calibrated by 199 permutations over 80 repeated trials for each pair of signal strength and sample size. Figure 9 supports the use of MS-SDMW as a kernel statistic for structured two-sample testing. At ∆ = 0, the rejection rate stays close to the nominal level. At ∆ = 0.04, power rises from 0.025 at eight metric spaces per group to 0.975 at 96 per group; at ∆ = 0.08, power reaches one by 32 samples per group. 8.7 Ablations The final experiments vary the multi-scale weights on MUTAG. We compare single orders, uniform averaging, inverse-order weights, and square-root inverse weights. Figure 10 is a useful boundary check. Multi-scale averaging is theoretically natural and stabilizes across orders, but it is not an oracle. On MUTAG, a single low order performs best. This does not contradict the population consistency statement in Theorem 5.22. That theorem concerns an asymptotic regime in which the scale range grows and the weights eventually move mass toward high orders. The MUTAG ablation is a finite-sample, finitegraph experiment in which high-order matrix laws are harder to estimate and local molecular patterns are highly informative. Thus inverse-order or low-order weights can be preferable in practice even though pure population consistency requires mass to escape to higher orders. 28
Distance-Matrix Wasserstein Statistics
Multi-scale weighting is helpful but dataset dependent
1.0
Nested-CV accuracy
0.9
0.8
0.7
0.6
0.5
le
g sin
2
n=
le
g sin
3
n=
le
g sin
4
n=
le
g sin
6
n=
le
g sin
8
n=
er
rm
ifo
un
inv
e ers
ord
inv qrt
e
ers
s
Figure 10: Multi-scale weight ablation on MUTAG. Low-order statistics are very strong on this small molecular dataset; multi-scale weights are useful defaults but should be tuned when labels are available.
A practical protocol is to use inverse-order weights as an unsupervised default and tune weights by nested validation in supervised tasks. 8.8 Empirical Takeaways The expanded experiments support a sharper conclusion than the preliminary study. DMW has a clear statistical role as a finite-subspace approximation to GW, and sliced multi-scale DMW has the expected computational behavior. It is particularly attractive for metric structural comparison and two-sample testing. At the same time, the graph classification table exposes an important limitation: metric-only DMW does not automatically recover the combinatorial information used by WL or graphlet kernels. This is not a failure of the theory; it identifies the regime where fused or label-aware DMW variants should be developed next.
9 Related Work DMW sits at the intersection of four lines of work: intrinsic optimal transport for metric spaces, random distance-matrix invariants, scalable relaxations of optimal transport, and graph substructure statistics. The closest mathematical relatives are GW lower bounds and reconstruction results based on distance distributions; the closest computational relatives are entropic and sliced GW methods. The distinction we emphasize throughout this section is where the relaxation is made. Entropic and sliced GW still operate on the original pointlevel alignment problem, either by regularizing it or by projecting it. DMW instead changes the object being transported: it compares laws of random finite metric subspaces. 29
Xu and Wu
Gromov–Wasserstein and structured optimal transport. This work builds on Gromov– Wasserstein distances and the metric approach to object matching (Gromov, 1999; Mémoli, 2011), computational optimal transport (Peyré and Cuturi, 2019), and optimal transport methods for structured data such as graphs (Vayer et al., 2019a). GW compares two metric measure spaces by optimizing over point-level couplings and measuring pairwise distance distortion. DMW keeps the same intrinsic viewpoint but replaces direct point-level alignment by optimal transport between laws of random finite distance matrices. The lower-bound theorem shows that this is a relaxation of GW, while the reverse approximation theorem shows that the relaxation becomes tight at the population level as the distance-matrix order grows. Entropic GW and regularized solvers. Entropic regularization is a central computational device in optimal transport, where Sinkhorn iterations turn large linear programs into matrix-scaling problems (Cuturi, 2013; Peyré and Cuturi, 2019). In the GW setting, entropically regularized objectives and related barycenter formulations provide practical algorithms for matching kernel or distance matrices (Peyré et al., 2016). These methods are important baselines in our experiments because they preserve the original GW modeling choice: they still optimize over a coupling between the points of the two input spaces, but make that optimization smoother and more numerically tractable. DMW addresses a different bottleneck. It removes the need to solve a point-level coupling for every object pair by sampling finite subspaces and transporting their distance-matrix laws. Thus entropic GW and DMW are complementary: entropic GW is a regularized solver for the alignment problem, whereas DMW is a statistical relaxation of the alignment problem. Distance matrix distributions, Gromov-weak topology, and reconstruction. Random distance matrices are a classical invariant of metric measure spaces. Gromov’s reconstruction principle states that the law of the infinite random distance array generated by i.i.d. samples determines a metric measure space up to measure-preserving isometry (Gromov, 1999). Related viewpoints appear in Vershik’s work on random metric spaces and universality (Vershik, 2004). In probability theory, finite-dimensional distance matrix distributions play a central role in the Gromov-weak topology and convergence in distribution of random metric measure spaces (Greven et al., 2009). Our full-hierarchy completeness proposition is an application of this existing reconstruction theory. The new element here is not reconstruction itself, but the use of Wasserstein distances between finite distance-matrix laws as scalable statistics with explicit GW approximation guarantees. Distance distributions and GW lower bounds. The order-2 case of DMW is closely related to distance distributions and inverse problems for metric measure spaces (Mémoli and Needham, 2022). It also connects to lower bounds for GW based on global distance distributions, shape distributions, and local distributions of distances. These lower bounds already show that comparing distributions of distances can provide computationally cheaper (n) GW proxies. Our higher-order construction differs in two ways. First, ρX records the joint law of all pairwise distances among n sampled points, so it captures triangles, four-point configurations, and higher-order finite metric patterns rather than only one-dimensional distance histograms. Second, our main theorem proves a reverse approximation inequality: 30
Distance-Matrix Wasserstein Statistics
the population DMW hierarchy not only lies below GW but converges back to GW with a gap controlled by empirical Wasserstein approximation of the underlying measures. Sliced, projected, and robust optimal transport. Sliced Wasserstein methods reduce high-dimensional OT to one-dimensional projections (Bonneel et al., 2015). Projection robust Wasserstein distances instead optimize over low-dimensional projections in order to retain discriminative structure while avoiding the full ambient dimension (Lin et al., 2020). Sliced Gromov–Wasserstein methods apply related projection ideas to GW-type structural comparison (Vayer et al., 2019b). Sliced DMW is different in where the projection enters: we first sample finite metric subspaces and map them to distance-matrix vectors, and only then slice the resulting distance-matrix laws. Thus the one-dimensional OT computations are performed on projected finite-subspace statistics, not directly on the original point clouds or on a GW objective. This distinction is important statistically: the sliced estimator targets SDMWn,p , a further relaxation of DMWn,p , and finite L controls only the Monte Carlo approximation of this sliced population statistic. Graph kernels and substructure statistics. For graph learning, DMW is related to graph kernels and substructure-counting methods. Shortest-path kernels compare distributions of graph distances (Borgwardt and Kriegel, 2005); graphlet kernels compare counts of small induced subgraphs (Shervashidze et al., 2009); and Weisfeiler–Lehman kernels compare iterative color-refinement features (Shervashidze et al., 2011). DMW can be interpreted as an optimal-transport comparison of random metric substructure distributions. Unlike graphlet or WL features, it is formulated for general metric measure spaces and comes with an explicit relation to GW. Unlike GW, it does not return a pointwise alignment, which is a limitation for matching tasks but an advantage for scalable structural comparison, kernels, and two-sample testing.
10 Discussion DMW is most useful when structural similarity is more important than recovering a pointwise alignment. This distinction is more than an implementation detail. In many learning tasks, including retrieval, clustering, kernel classification, and two-sample testing, the final object used by the algorithm is a scalar dissimilarity or a kernel value, not the coupling matrix itself. Standard GW nevertheless spends most of its computational effort searching for a global point-level alignment. DMW asks what can be gained by giving up that alignment as an output and comparing the distribution of finite metric patterns instead. The results in this paper show that this change of viewpoint is mathematically controlled: at the population level, finite-subspace laws form a hierarchy that converges back to GW, while at the algorithmic level the sampling parameters expose an explicit statistical–computational tradeoff. The framework also leaves several theoretical and computational questions open. Fixedorder DMW is not complete, full empirical DMW can suffer from high-dimensional Wasserstein estimation, and sliced DMW is a further relaxation introduced for scalability. The intrinsic-dimension bounds and the sliced W1 estimates explain why the method can work beyond the ambient dimension of the distance-matrix vector, but they are not the end of the story. One natural direction is to replace purely random slicing by more adaptive pro31
Xu and Wu
jections, such as projection-robust or learned low-dimensional summaries, while preserving the lower-bound and consistency structure. Another is to understand optimal multi-scale weighting beyond the simple oracle bound in this paper. In finite samples, the best weights depend on the approximation error at each order, the estimation variance of the corresponding matrix law, and the downstream task. A fully adaptive theory of scale selection for distance-matrix laws remains open. The experiments also clarify an important modeling limitation. The DMW variant studied here is deliberately metric: graphs are represented through shortest-path distances and a node measure, without using node labels or attributes except through the chosen metric. This is the right abstraction for intrinsic structural comparison, but it does not capture all information used by modern graph kernels and graph neural networks. The gap to Weisfeiler–Lehman and graphlet kernels on some TU datasets is therefore informative rather than surprising. It points to a natural next step: a fused DMW statistic that combines distance-matrix laws with feature or label laws, in the same spirit that fused GW augments structural distortion with feature matching. Such a construction would preserve the finite-subspace perspective while allowing the sampled subspaces to carry both geometry and attributes.
Acknowledgments and Disclosure of Funding This work is supported by the Zhongguancun Academy, Grant No. C20250201.
Appendix A. Mathematical Background Tools For completeness and to make the proof dependencies explicit, we collect the standard results used throughout the appendix. None of the statements in this section is new; they are included only to clarify where concentration, projection, and reconstruction inputs enter the argument. Lemma A.1 (McDiarmid’s bounded-difference inequality (McDiarmid, 1989)) Let Z1 , . . . , Zm be independent random variables and let F = F (Z1 , . . . , Zm ). If changing only Zi can change F by at most ci , then for every t > 0, 2t2 Pr{|F − EF | ≥ t} ≤ 2 exp − Pm 2 . i=1 ci Lemma A.2 (Hoeffding’s inequality (Hoeffding, 1963)) Let Z1 , . . . , ZL be independent random variables with 0 ≤ Zℓ ≤ M . Then, for every t > 0, ( ) L 1X 2Lt2 Pr Zℓ − EZ1 ≤ −t ≤ exp − 2 . L M ℓ=1
Lemma A.3 (One-dimensional empirical Wasserstein (Dvoretzky et al., 1956; Massart, 1990)) The DKW inequality implies the following bound. Let α be supported on an interval of length at most R and let α bK be its empirical measure from K i.i.d. samples. Then EW1 (b αK , α) ≤ CRK −1/2 32
Distance-Matrix Wasserstein Statistics
for a universal constant C. Lemma A.4 (Empirical Wasserstein rates (Fournier and Guillin, 2015; Weed and Bach, 2019)) Let α be a compactly supported probability measure on a metric space (Z, d), and assume that its support satisfies a polynomial covering bound N (supp α, d, ε) ≤ Aε−s ,
0 < ε ≤ 1.
Then standard compact-support empirical Wasserstein bounds, obtained for example from the entropy methods in Fournier and Guillin (2015) and Weed and Bach (2019), imply −1/(2p) , K EWp (b αK , α) ≤ C K −1/(2p) (log K)1/p , −1/s K ,
s < 2p, s = 2p, s > 2p,
with C depending on the support diameter, A, and (s, p). Lemma A.5 (Cramer–Wold theorem (Cramér and Wold, 1936; Billingsley, 1999)) Two Borel probability measures on Rd are equal if all of their one-dimensional linear projections have the same law. Moreover, if αm , α ∈ P(Rd ) and (Pθ )# αm ⇒ (Pθ )# α for every θ ∈ Rd , then αm ⇒ α. Lemma A.6 (Gromov reconstruction principle (Gromov, 1999; Greven et al., 2009)) For compact metric measure spaces, the joint laws of all finite random distance matrices determine the space up to measure-preserving isometry. Lemma A.7 (Schoenberg criterion (Schoenberg, 1938; Berg et al., 1984)) If a symmetric function d is conditionally negative definite, then exp(−λd) is positive definite for every λ > 0.
Appendix B. Proofs of Auxiliary Results B.1 Proof of Lemma B.1 Lemma B.1 (Triangle inequality for GW) Let X = (X, dX , µ), Y = (Y, dY , ν), and Z = (Z, dZ , ζ) be compact metric measure spaces. Then GWp (X , Z) ≤ GWp (X , Y) + GWp (Y, Z). Proof Let πXY ∈ Π(µ, ν) and πY Z ∈ Π(ν, ζ). By the gluing lemma, there exists a probability measure λ on X × Y × Z whose (X, Y ) marginal is πXY and whose (Y, Z) marginal is πY Z . Let πXZ be the (X, Z) marginal of λ. Then πXZ ∈ Π(µ, ζ), so it is admissible for GWp (X , Z). Let (x, y, z) and (x′ , y ′ , z ′ ) be two independent samples from λ. The triangle inequality in R gives |dX (x, x′ ) − dZ (z, z ′ )| ≤ |dX (x, x′ ) − dY (y, y ′ )| + |dY (y, y ′ ) − dZ (z, z ′ )|. 33
Xu and Wu
Taking the Lp (λ ⊗ λ) norm and applying Minkowski’s inequality yields 1/p 1/p 1/p E|dX (x, x′ ) − dZ (z, z ′ )|p ≤ E|dX (x, x′ ) − dY (y, y ′ )|p + E|dY (y, y ′ ) − dZ (z, z ′ )|p . The first term on the right depends only on πXY , and the second depends only on πY Z . Therefore GWp (X , Z) ≤ Dp (πXY ) + Dp (πY Z ), where Dp denotes the corresponding GW distortion functional. Taking the infimum over πXY and πY Z gives the claim.
B.2 Proof of Lemma B.2 Lemma B.2 (Same-space stability) Let (X, d) be a compact metric space. For probability measures µ, µ′ on X, GWp (X, d, µ), (X, d, µ′ ) ≤ 2Wp (µ, µ′ ). Proof Let γ ∈ Π(µ, µ′ ). Use γ as an admissible GW coupling between (X, d, µ) and (X, d, µ′ ). For (x, y), (x′ , y ′ ) ∈ X × X, the triangle inequality gives |d(x, x′ ) − d(y, y ′ )| ≤ d(x, y) + d(x′ , y ′ ) because d(x, x′ ) ≤ d(x, y) + d(y, y ′ ) + d(y ′ , x′ ) and symmetrically with (x, x′ ) and (y, y ′ ) interchanged. Therefore, by Minkowski’s inequality in Lp (γ ⊗ γ), Z Z 1/p Z 1/p ′ ′ p ′ ′ p |d(x, x ) − d(y, y )| dγ(x, y) dγ(x , y ) ≤2 d(x, y) dγ(x, y) . Minimizing over γ proves the result.
B.3 Proof of Lemma B.3 (n)
(n)
Lemma B.3 (Lifting lemma) Let Γ ∈ Π(ρX , ρY ). There exists a probability measure Λ on X n × Y n whose marginals are µ⊗n and ν ⊗n and such that (n)
(n)
(TX , TY )# Λ = Γ. Proof Since X n and Y n are compact metric spaces, they are standard Borel spaces. Hence regular conditional distributions exist. Let {αa : a ∈ RNn },
{βb : b ∈ RNn } (n)
(n)
be regular conditional distributions of µ⊗n and ν ⊗n given the maps TX and TY , respectively. Thus, for every Borel set A ⊆ X n , Z (n) ⊗n µ (A) = αa (A) dρX (a), RNn
34
Distance-Matrix Wasserstein Statistics
(n)
(n)
and αa is concentrated on the fiber (TX )−1 ({a}) for ρX -almost every a. Analogous state(n) ments hold for βb and ρY . Define a probability measure Λ on X n × Y n by Z αa (A) βb (B) dΓ(a, b) Λ(A × B) = RNn ×RNn
for Borel rectangles A ⊆ X n and B ⊆ Y n , and extend it to the product σ-algebra. Its first marginal is, for Borel A ⊆ X n , Z Z (n) n Λ(A × Y ) = αa (A) dΓ(a, b) = αa (A) dρX (a) = µ⊗n (A), (n)
where we used that the first marginal of Γ is ρX . The second marginal is similarly ν ⊗n . It remains to check the pushforward identity. For Borel sets C, D ⊆ RNn , (n) (n) (n) (n) (TX , TY )# Λ(C × D) = Λ (TX )−1 (C) × (TY )−1 (D) . Using the fiber concentration of αa and βb , the last display equals Z χC (a)χD (b) dΓ(a, b) = Γ(C × D). (n)
(n)
Since rectangles generate the Borel σ-algebra on RNn × RNn , we obtain (TX , TY )# Λ = Γ.
B.4 Proof of Lemma B.4 (n)
(n)
Lemma B.4 (Empirical GW controlled by DMW) Let Γ⋆ ∈ Π(ρX , ρY ) be optimal, and let Λ⋆ be any lift of Γ⋆ given by Lemma B.3. For a sample ((x1 , . . . , xn ), (y1 , . . . , yn )) ∼ Λ⋆ , define n
µ bn =
n
1X δxi , n
νbn =
i=1
1X δyi . n i=1
Then EΛ⋆ GWp (Xbn , Ybn ) ≤
n−1 n
1/p DMWn,p (X , Y).
Proof Fix a pair of tuples (x1 , . . . , xn , y1 , . . . , yn ). Consider the empirical coupling n
πn =
1X δ(xi ,yi ) ∈ Π(b µn , νbn ). n i=1
35
Xu and Wu
Using πn as an admissible coupling in the definition of GW gives n X
1 GWpp (Xbn , Ybn ) ≤ 2 n
|dX (xi , xj ) − dY (yi , yj )|p .
i,j=1
The diagonal terms vanish, and the two off-diagonal orientations contribute equally. Hence GWpp (Xbn , Ybn ) ≤
2 n2
X
|dX (xi , xj ) − dY (yi , yj )|p .
1≤i<j≤n
Since Nn = n(n − 1)/2, the factor 2/n2 is equal to ((n − 1)/n)(1/Nn ). Hence n−1 1 GWpp (Xbn , Ybn ) ≤ n Nn
X
|dX (xi , xj ) − dY (yi , yj )|p .
1≤i<j≤n (n)
(n)
Taking expectation under Λ⋆ and using (TX , TY )# Λ⋆ = Γ⋆ yields Z n−1 n−1 p b b ⋆ EΛ GWp (Xn , Yn ) ≤ ∥a − b∥pp,avg dΓ⋆ (a, b) = DMWpn,p (X , Y). n n Finally, Jensen’s inequality gives EΛ⋆ GWp (Xbn , Ybn ) ≤
EΛ⋆ GWpp (Xbn , Ybn )
1/p
≤
n−1 n
1/p DMWn,p (X , Y).
Appendix C. Proofs Deferred from the Main Text C.1 Proof of Proposition 3.2 Proof By definition, (n)
(k)
Pn→k ◦ TX (x1 , . . . , xn ) = TX (x1 , . . . , xk ). Since the first k coordinates of µ⊗n have law µ⊗k , the pushforward identity follows. The same argument applies to any k-element subset S ⊂ {1, . . . , n}: if PS : RNn → RNk keeps the distances among the indices in S, then (n)
(k)
(n)
(PS )# ρX = ρX ,
(k)
(PS )# ρY = ρY .
(n)
(n)
Let Γ be an optimal coupling of ρX and ρY . For every S, the measure (PS , PS )# Γ is (k) (k) a coupling of ρX and ρY , hence Z DMWpk,p (X , Y) ≤ ∥PS a − PS b∥pp,avg,k dΓ(a, b). 36
Distance-Matrix Wasserstein Statistics
Average this inequality uniformly over all k-element subsets S. For fixed a, b ∈ RNn , each n pair 1 ≤ i < j ≤ n appears in exactly n−2 k−2 of the k subsets. Therefore 1 X 1 ∥PS a − PS b∥pp,avg,k = n Nn k |S|=k
X
|aij − bij |p = ∥a − b∥pp,avg,n .
1≤i<j≤n
Integrating over Γ gives DMWpk,p (X , Y) ≤
Z
∥a − b∥pp,avg,n dΓ(a, b) = DMWpn,p (X , Y).
Taking pth roots proves the monotonicity claim.
C.2 Proof of Proposition 3.4 Proof Fix π ∈ Π(µ, ν) and set (n)
(n)
Γπ = (TX , TY )# π ⊗n . The first marginal of π ⊗n on X n is µ⊗n , and the second marginal on Y n is ν ⊗n . Therefore the first marginal of Γπ is (n) (n) (TX )# µ⊗n = ρX , and the second marginal is (n)
(n)
(TY )# ν ⊗n = ρY . (n)
(n)
Thus Γπ ∈ Π(ρX , ρY ). Since this holds for every π ∈ Π(µ, ν), the inclusion follows.
C.3 Proof of Theorem 4.1 Proof Fix π ∈ Π(µ, ν) and draw (xi , yi )ni=1 independently from π. Let Λπ = π ⊗n be the resulting probability measure on (X × Y )n , identified with X n × Y n . The pushforward (n)
(n)
Γπ = (TX , TY )# Λπ (n)
(n)
(n)
is a coupling of ρX and ρY , because the first marginal is (TX )# µ⊗n and the second (n) marginal is (TY )# ν ⊗n . Hence Z (n) (n) DMWpn,p (X , Y) ≤ ∥TX (x1 , . . . , xn ) − TY (y1 , . . . , yn )∥pp,avg dΛπ . By the definition of ∥ · ∥p,avg and Fubini’s theorem, the right-hand side is X Z 1 |dX (xi , xj ) − dY (yi , yj )|p dΛπ . Nn 1≤i<j≤n
37
Xu and Wu
Since the pairs (xi , yi ) are independent with law π, each summand is equal to ZZ |dX (x, x′ ) − dY (y, y ′ )|p dπ(x, y) dπ(x′ , y ′ ). There are Nn identical summands, so ZZ p DMWn,p (X , Y) ≤ |dX (x, x′ ) − dY (y, y ′ )|p dπ(x, y) dπ(x′ , y ′ ). Taking the infimum over π ∈ Π(µ, ν) gives DMWpn,p (X , Y) ≤ GWpp (X , Y). The claim follows by taking pth roots.
C.4 Proof of Theorem 4.2 Proof Let Γ⋆ and Λ⋆ be as in Lemma B.4. Under Λ⋆ , the X n marginal is µ⊗n and the Y n marginal is ν ⊗n . Thus µ bn and νbn are ordinary empirical measures built from n i.i.d. samples from µ and ν, respectively, although they may be dependent across the two spaces. For every realization of the two empirical spaces, Lemma B.1 gives GWp (X , Y) ≤ GWp (X , Xbn ) + GWp (Xbn , Ybn ) + GWp (Ybn , Y). Taking expectation with respect to Λ⋆ gives GWp (X , Y) ≤ EΛ⋆ GWp (X , Xbn ) + EΛ⋆ GWp (Xbn , Ybn ) + EΛ⋆ GWp (Ybn , Y). Lemma B.2 bounds the first and third terms by 2EWp,X (µ, µ bn ),
2EWp,Y (ν, νbn ),
respectively. The middle term is bounded by ((n − 1)/n)1/p DMWn,p (X , Y) by Lemma B.4. Since ((n − 1)/n)1/p ≤ 1, combining the three bounds proves the stated inequality.
C.5 Proof of Corollary 4.3 Proof The lower bound follows from Theorem 4.1. The upper bound follows by rearranging Theorem 4.2, using ((n − 1)/n)1/p ≤ 1.
C.6 Proof of Corollary 4.4 Proof On a compact metric space, empirical measures converge almost surely to the population measure in Wp , and Wp (µ, µ bn ) is bounded by diam(X). Dominated convergence gives EWp (µ, µ bn ) → 0, and similarly for ν. Corollary 4.3 then implies 0 ≤ GWp (X , Y) − DMWn,p (X , Y) → 0.
38
Distance-Matrix Wasserstein Statistics
C.7 Proof of Corollary 4.5 Proof Let X = {x1 , . . . , xm } and write pi = µ({xi }) and pbi = µ bn ({xi }). Since all unmatched mass can be transported at cost at most RX , m
p p Wp,X (µ, µ bn ) ≤ RX TV(µ, µ bn ),
TV(µ, µ bn ) =
1X |b pi − pi |. 2 i=1
By Cauchy–Schwarz and the multinomial variance identity,
E
m X i=1
v r u m r P u X m(1 − i p2i ) m 2 t |b pi − pi | ≤ m E(b p i − pi ) = ≤ . n n i=1
√ √ Thus E TV(µ, µ bn ) ≤ m/(2 n). Jensen’s inequality applied to the concave map t 7→ t1/p gives the claimed bound for X. The proof for Y is identical, and the GW–DMW bound follows from Corollary 4.3.
C.8 Proof of Theorem 5.1 Proof For any four probability measures α, α′ , β, β ′ on a common metric space, the triangle inequality implies Wp (α, β) ≤ Wp (α, α′ ) + Wp (α′ , β ′ ) + Wp (β ′ , β) and, by the same argument with primed and unprimed measures exchanged, Wp (α′ , β ′ ) ≤ Wp (α′ , α) + Wp (α, β) + Wp (β, β ′ ). Therefore |Wp (α, β) − Wp (α′ , β ′ )| ≤ Wp (α, α′ ) + Wp (β, β ′ ). Apply this deterministic inequality with (n)
α = ρbX,K ,
(n)
β = ρbY,K ,
(n)
α′ = ρX ,
(n)
β ′ = ρY .
Since |EZ − c| ≤ E|Z − c| for any integrable random variable Z and constant c, taking expectations yields (n)
(n)
(n)
(n)
\ n,p,K − DMWn,p (X , Y) ≤ EWp (b EDMW ρX,K , ρX ) + EWp (b ρY,K , ρY ), which is the desired bound.
39
Xu and Wu
C.9 Proof of Theorem 5.2 Proof Let (n)
(n)
F (DX,1 , . . . , DX,K , DY,1 , . . . , DY,K ) = Wp (b ρX,K , ρbY,K ). (n)
(n)
The supports of ρX and ρY have diameter at most R in ∥ · ∥p,avg , because all distance entries lie in intervals of length at most R. If one X-sample is replaced, the empirical (n) measure ρbX,K changes by replacing one atom of mass 1/K. Coupling all unchanged atoms to themselves and the replaced atom to the new atom gives (n)
(n)′
Wp (b ρX,K , ρbX,K ) ≤ RK −1/p . The same bound holds when one Y -sample is replaced. By the reverse triangle inequality for Wasserstein distance, changing any one of the 2K input variables changes F by at most RK −1/p . By Lemma A.1, for every t > 0, 2t2 t2 Pr{|F − EF | ≥ t} ≤ 2 exp − = 2 exp − 2 1−2/p . 2KR2 K −2/p R K With probability at least 1 − δ, r |F − EF | ≤ RK
1/2−1/p
2 log . δ
Finally, Theorem 5.1 gives |EF − DMWn,p (X , Y)| ≤ ηX (n, K) + ηY (n, K), and the result follows by the triangle inequality on the real line.
C.10 Proof of Theorem 5.4 Proof Let (n)
(n)
G(DX,1 , . . . , DX,K , DY,1 , . . . , DY,K ) = Wpp (b ρX,K , ρbY,K ). The cost ∥a − b∥pp,avg is bounded by Rp on the product of the two distance-matrix supports. (n) If one atom of ρbX,K is replaced, then only mass 1/K must be reassigned relative to any feasible transport plan. Since the cost is bounded by Rp , the optimal transport cost changes by at most Rp /K. The same argument applies when one Y atom is replaced. Thus G has bounded differences with constants Rp /K over 2K independent inputs. By Lemma A.1, 2t2 Kt2 Pr{|G − EG| ≥ t} ≤ 2 exp − = 2 exp − 2p . 2K(Rp /K)2 R Solving for t at confidence level δ proves the claim.
40
Distance-Matrix Wasserstein Statistics
C.11 Proof of Theorem 5.6 Proof By the triangle inequality on the real line, \ n,p,K − GWp (X , Y) ≤ EDMW \ n,p,K − DMWn,p (X , Y) +|DMWn,p (X , Y) − GWp (X , Y)| . EDMW The first term is bounded by Theorem 5.1. The second term is bounded by Corollary 4.3.
C.12 Proof of Corollary 5.7 Proof By the triangle inequality, \ n,p,K − GWp (X , Y)| ≤ |DMW \ n,p,K − DMWn,p (X , Y)| + |DMWn,p (X , Y) − GWp (X , Y)|. |DMW The first term is bounded with probability at least 1 − δ by Theorem 5.2. The second term is bounded deterministically by Corollary 4.3.
C.13 Proof of Proposition 5.8 Proof Under the null hypothesis, the joint law of the pooled sample is invariant under all permutations of the group labels that preserve the group sizes. Conditional on the unordered pooled sample, the observed labeling is therefore uniformly distributed over the permutation orbit. The exact permutation p-value is the fraction of labelings whose statistic is at least as large as the observed statistic, with any ties handled conservatively or by randomized tie-breaking. Hence, conditional on the pooled data, this p-value is stochastically no smaller than a uniform random variable on [0, 1]. Therefore the conditional rejection probability is at most α, and integrating over the pooled sample gives unconditional type-I error at most α.
C.14 Proof of Corollary 5.9 Proof Apply Lemma A.4, in its Euclidean compact-support form (Fournier and Guillin, (n) (n) 2015), to the compactly supported laws ρX and ρY on RNn .
C.15 Proof of Proposition 5.14 (n)
Proof For each θ ∈ Θn,p , the map Pθ is 1-Lipschitz. Therefore, for any coupling Γ of ρX (n) and ρY , Z Z (n) (n) Wpp (Pθ )# ρX , (Pθ )# ρY ≤ |Pθ (a) − Pθ (b)|p dΓ(a, b) ≤ ∥a − b∥pp,avg dΓ(a, b). Taking the infimum over Γ gives (n)
(n)
Wp (Pθ )# ρX , (Pθ )# ρY
41
≤ DMWn,p (X , Y).
Xu and Wu
Integrating the pth powers over θ ∼ σ yields SDMWn,p ≤ DMWn,p . The second inequality is Theorem 4.1.
C.16 Proof of Theorem 5.18 Proof For each θ ∈ Θn,p , the projected supports have diameter at most R because Pθ is 1-Lipschitz and the distance-matrix supports have diameter at most R under ∥ · ∥p,avg . Hence the random variable (n) (n) Zθ = Wpp (Pθ )# ρbX,K , (Pθ )# ρbY,K p
\ n,p,K,L is the takes values in [0, Rp ]. Conditional on the empirical distance matrices, SDMW p average of L i.i.d. copies of Zθ , with mean SDMWn,p,K . Therefore v ! u L u X p Rp 1 p t \ n,p,K,L − SDMW Zθℓ ≤ √ , E SDMW Var n,p,K ≤ L 2 L ℓ=1 using the variance bound for a random variable supported on an interval of length Rp .
C.17 Proof of Corollary 5.19 Proof Let (n)
(n)
Zℓ = Wpp (Pθℓ )# ρX , (Pθℓ )# ρY
.
Then 0 ≤ Zℓ ≤ Rp and EZℓ = SDMWpn,p (X , Y) ≥ εp . By Lemma A.2, ) ( L 2Lγ 2 ε2p 1X p Zℓ < EZℓ − γε ≤ exp − . Pr L R2p ℓ=1
The stated sufficient condition on L makes this probability at most δ. Since EZℓ − γεp ≥ (1 − γ)εp , the claim follows.
C.18 Proof of Proposition 5.20 Proof For each θ, the reverse triangle inequality for Wasserstein distance gives (n)
(n)
(n)
(n)
(n)
(n)
(n)
(n)
Wp ((Pθ )# ρbX,K , (Pθ )# ρbY,K ) − Wp ((Pθ )# ρX , (Pθ )# ρY ) ≤ Wp ((Pθ )# ρbX,K , (Pθ )# ρX ) + Wp ((Pθ )# ρbY,K , (Pθ )# ρY ). Since Pθ is 1-Lipschitz, the two terms on the right are bounded by the corresponding unprojected Wasserstein errors, (n)
(n)
Wp (b ρX,K , ρX )
and 42
(n)
(n)
Wp (b ρY,K , ρY ).
Distance-Matrix Wasserstein Statistics
All relevant Wasserstein distances are bounded by R, so |ap − bp | ≤ pRp−1 |a − b| for a, b ∈ [0, R]. Integrating over θ and taking expectation over the empirical distance matrices gives the first claim. The one-dimensional K −1/2 bound for W1 follows from Lemma A.3. This is the statistical advantage motivating sliced DMW: each projected problem is one-dimensional even when Nn is large.
C.19 Proof of Theorem 5.21 Proof For p = 1, the sliced statistic is linear in the directional Wasserstein distances: (n)
(n)
SDMWn,1 (X , Y) = Eθ∼σ W1 (Pθ )# ρX , (Pθ )# ρY
.
Define the empirical directional average (n) (n) SDMWn,1,K = Eθ∼σ W1 (Pθ )# ρbX,K , (Pθ )# ρbY,K .
By the triangle inequality, \ n,1,K,L −SDMWn,1 | ≤ E|SDMW \ n,1,K,L −SDMWn,1,K |+E|SDMWn,1,K −SDMWn,1 |. E|SDMW √ The first term is bounded by R/(2 L) by Theorem 5.18 with p = 1. For the second term, use the reverse triangle inequality for W1 and integrate over θ: (n)
(n)
|SDMWn,1,K − SDMWn,1 | ≤ Eθ∼σ W1 (Pθ )# ρbX,K , (Pθ )# ρX (n)
(n)
+Eθ∼σ W1 (Pθ )# ρbY,K , (Pθ )# ρY
.
Taking expectation over the empirical distance matrices and applying the assumed onedimensional rates gives 2C1 RK −1/2 .
C.20 Proof of Proposition 7.1 Proof Let X = {1, 2, 3, 4} with uniform measure, and define dX by declaring the pairs (1, 3) and (2, 4) to have distance 2, all other distinct pairs to have distance 1, and dX (i, i) = 0. This is the shortest-path metric on the four-cycle. Let Y = {1, 2, 3, 4} with uniform measure, and define dY by declaring the pairs (1, 2) and (1, 3) to have distance 2, all other distinct pairs to have distance 1, and dY (i, i) = 0. The triangle inequality holds: any distance equal to 2 is realized as the sum of two unit distances through the remaining points, and all other nontrivial triangle inequalities are immediate. In both spaces, among the six unordered pairs of distinct points, exactly four pairs have distance 1 and exactly two pairs have distance 2. Since the measures are uniform, the law of d(x1 , x2 ) for two i.i.d. samples is the same in X and Y : both have mass 1/4 at 0, mass (2) (2) 1/2 at 1, and mass 1/4 at 2. Hence ρX = ρY . 43
Xu and Wu
The spaces are not isometric. Indeed, in X the graph whose edges connect pairs at distance 2 is a matching, so every point has degree 1 in this graph. In Y , the distance-2 graph has edges (1, 2) and (1, 3), so one point has degree 2 and one point has degree 0. This degree pattern is invariant under isometry, proving non-isometry.
C.21 Proof of Proposition 7.2 (n)
Proof The family {ρX }n≥2 is precisely the family of finite-dimensional distributions of the infinite random distance matrix i.i.d. dX (xi , xj ) 1≤i<j<∞ , x1 , x2 , . . . ∼ µ. Equality of all finite-dimensional distributions for X and Y implies equality in law of the corresponding infinite random distance arrays. By Lemma A.6, the law of this infinite distance array determines the metric measure space up to measure-preserving isometry. Therefore X and Y are measure-preserving isometric.
C.22 Proof of Theorem 5.10 Proof We give the argument for X; the proof for Y is identical. Equip X n with the normalized product metric n
dX,n,p ((xi ), (x′i )) =
1X dX (xi , x′i )p n
!1/p .
i=1
For two tuples x = (xi ) and x′ = (x′i ), the reverse triangle inequality gives |dX (xi , xj ) − dX (x′i , x′j )| ≤ dX (xi , x′i ) + dX (xj , x′j ). Using (a + b)p ≤ 2p−1 (ap + bp ) and counting how often each index appears among unordered pairs, (n)
(n)
∥TX (x) − TX (x′ )∥pp,avg ≤
2p−1 X {dX (xi , x′i )p + dX (xj , x′j )p } = 2p dX,n,p (x, x′ )p . Nn i<j
(n)
Thus TX is 2-Lipschitz from (X n , dX,n,p ) into (RNn , ∥ · ∥p,avg ). If {z1 , . . . , zM } is an ε-net of supp µ in dX , then the product grid is an ε-net of (supp µ)n in dX,n,p . Therefore N ((supp µ)n , dX,n,p , ε) ≤ N (supp µ, dX , ε)n ≤ AnX ε−nsX . (n)
Since TX is 2-Lipschitz, (n) N (supp ρX , ∥ · ∥p,avg , ε) ≤ N ((supp µ)n , dX,n,p , ε/2) ≤ AnX
44
nsX 2 . ε
Distance-Matrix Wasserstein Statistics
This proves the stated covering estimate. Finally, empirical Wasserstein bounds on compact metric spaces with polynomial covering dimension s give (n)
(n)
EWp (b ρX,K , ρX ) ≤ CX (n, p)ψnsX ,p (K), with the displayed cases for ψs,p ; this is exactly Lemma A.4. Applying the same argument to Y completes the proof.
C.23 Proof of Corollary 5.12 Proof Compact Riemannian manifolds with the stated bounded-geometry assumptions have polynomial covering numbers at their intrinsic dimension: there exists AX such that N (X, dX , ε) ≤ AX ε−sX for 0 < ε ≤ 1, and similarly for Y . This follows from the standard volume comparison and normal-coordinate covering argument: the lower injectivity-radius bound gives uniformly controlled coordinate balls, the curvature and diameter bounds control the number of such balls needed at scale ε, and compactness handles large scales. The result is then an immediate application of Theorem 5.10.
C.24 Proof of Proposition 5.15 (n)
(n)
Proof For fixed n, ρX and ρY are Borel probability measures on the finite-dimensional vector space RNn . If all one-dimensional projections agree, then for every θ ∈ RNn and every t ∈ R, Z Z (n)
eit⟨θ,a⟩ dρX (a) =
(n)
(n)
eit⟨θ,a⟩ dρY (a). (n)
Equivalently, the characteristic functions of ρX and ρY agree at every frequency ξ = tθ ∈ (n) (n) RNn . By the Cramer–Wold theorem collected in Lemma A.5, ρX = ρY . If this holds for every n ≥ 2, then the full distance-matrix hierarchies are equal, and Proposition 7.2 implies that X and Y are measure-preserving isometric.
C.25 Proof of Corollary 5.16 Proof Fix n. Since SDMWn,p = 0, the nonnegative function (n) (n) f (θ) = Wp (Pθ )# ρX , (Pθ )# ρY vanishes for σ-almost every θ. The map f is continuous on Θn,p : if the matrix laws are supported in a ball of radius Rn , then |f (θ) − f (θ′ )| ≤ 2Rn ∥θ − θ′ ∥q,dual , 45
Xu and Wu
because changing the projection direction moves each projected random variable by at most Rn ∥θ − θ′ ∥q,dual . Since σ has full support, a continuous nonnegative function that is zero σ-almost surely must vanish everywhere on Θn,p . For every t ∈ R, the map Z Z (n) (n) θ 7→ eit⟨θ,a⟩ dρX (a) − eit⟨θ,a⟩ dρY (a) therefore vanishes for every θ ∈ Θn,p . Since every nonzero direction in RNn is a positive scalar multiple of a vector in Θn,p , equality of projected laws holds in every linear direction. (n) (n) By Proposition 5.15, ρX = ρY . Since this holds for every n ≥ 2, Proposition 7.2 implies that X and Y are measure-preserving isometric.
C.26 Proof of Theorem 5.25 (n)
(n)
Proof Fix n and write αm = ρXm and α = ρX . Uniform boundedness of the diameters implies that all these laws are supported in a common compact ball Bn ⊂ RNn . For θ ∈ Θn,p , define fm (θ) = Wp (Pθ )# αm , (Pθ )# α . The assumption SDMWn,p (Xm , X ) → 0 says that Z fm (θ)p dσ(θ) → 0. Θn,p
The family {fm }m is uniformly Lipschitz. Indeed, if Bn is contained in a ball of radius Rn for ∥ · ∥p,avg , then |fm (θ) − fm (θ′ )| ≤ 2Rn ∥θ − θ′ ∥q,dual
for all m, θ, θ′ ∈ Θn,p .
If fm (θ0 ) failed to converge to zero at some θ0 , this uniform Lipschitz bound would give an open neighborhood on which fm stays bounded away from zero along a subsequence. Since σ has full support, that neighborhood has positive σ-mass, contradicting the convergence of the integral. Hence fm (θ) → 0 for every θ ∈ Θn,p . By scaling, the same convergence of one-dimensional projected laws holds for every linear direction in RNn . The Lemma A.5 then implies weak convergence (n)
(n)
ρXm ⇒ ρX . Because the supports are uniformly bounded, weak convergence is equivalent to convergence in Wasserstein distance on RNn for fixed n. Hence every finite distance-matrix marginal converges. The Gromov-weak topology is precisely the topology induced by convergence of all finite distance-matrix distributions. Therefore Xm → X in the Gromov-weak topology. Under the additional precompactness and topology-agreement assumption stated in the theorem, this convergence implies GWp (Xm , X ) → 0.
46
Distance-Matrix Wasserstein Statistics
C.27 Proof of Proposition 5.27 Proof We use Lemma A.7 (Schoenberg, 1938; Berg et al., 1984). For probability measures α, β on the real line with finite first moment, Z |Fα (t) − Fβ (t)| dt, W1 (α, β) = R
where Fα and Fβ are their distribution functions. Thus W1 is the L1 distance between 1 CDF P embeddings. The L metric is conditionally negative definite: for coefficients ci with i ci = 0, Z X X ci cj ∥fi − fj ∥L1 = ci cj |fi (t) − fj (t)| dt ≤ 0, i,j
i,j
because the absolute-value metric on R is conditionally negative definite and the integral preserves conditional negative definiteness. For fixed θ, the map (n)
(n)
(X , Y) 7→ W1 ((Pθ )# ρX , (Pθ )# ρY ) is therefore conditionally negative definite. This argument does not require Pθ to be 1Lipschitz; the compactness of the distance-matrix supports and finite first moment of the direction law ensure the integral is finite. Integrating over θ preserves conditional negative definiteness, so SDMWn,1 is conditionally negative definite. Nonnegative sums over n preserve conditional negative definiteness, hence MSDMW1 is conditionally negative definite. Lemma A.7 gives the positive definiteness of the exponential kernels.
C.28 Proof of Theorem 5.22 Proof For each n, Proposition 5.14 gives SDMWn,p (X , Y) ≤ DMWn,p (X , Y) ≤ GWp (X , Y). Multiplying by nonnegative weights αn and summing over n ∈ N proves MSDMWp ≤ MDMWp ≤ GWp , P using n∈N αn = 1. For the error bound, apply Corollary 4.3 at each order: 0 ≤ GWp (X , Y) − DMWn,p (X , Y) ≤ 2EWp,X (µ, µ bn ) + 2EWp,Y (ν, νbn ). P Since GWp − MDMWp = n∈N αn (GWp − DMWn,p ), the displayed multi-scale error bound follows. Finally, assume that αm,n is a sequence of probability weights on {2, 3, . . .} whose mass escapes to high orders. By Corollary 4.4, DMWn,p (X , Y) → GWp (X , Y). The quantities 47
Xu and Wu
GWp − DMWn,p are nonnegative and bounded by GWp . For any ε > 0, choose N such that GWp − DMWn,p ≤ ε for all n ≥ N . Then X X GWp − MDMWp,m ≤ GWp αm,n + ε αm,n . n<N
n≥N
Taking m → ∞ and then ε ↓ 0 proves MDMWp,m → GWp .
C.29 Proof of Proposition 5.24 Proof The first claim is immediate because a linear function over the probability simplex attains its minimum at an extreme point, and the minimizing extreme points are exactly the scales with minimal Bn . For the entropy-regularized problem, form the Lagrangian ! X X X L(α, λ) = αn Bn + τ αn − 1 . αn log αn + λ n∈N
n∈N
n∈N
The first-order condition for an interior optimum is Bn + τ (1 + log αn ) + λ = 0. P Thus αn = C exp(−Bn /τ ) for a normalizing constant C. Enforcing n αn = 1 gives the displayed softmin weights. Strict convexity of the entropy-regularized objective on the simplex gives uniqueness.
References Christian Berg, Jens Peter Reus Christensen, and Paul Ressel. Harmonic Analysis on Semigroups: Theory of Positive Definite and Related Functions. Springer, 1984. doi: 10.1007/978-1-4612-1128-0. Patrick Billingsley. Convergence of Probability Measures. Wiley, 2 edition, 1999. ISBN 9780471197454. doi: 10.1002/9780470316962. Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister. Sliced and radon wasserstein barycenters of measures. Journal of Mathematical Imaging and Vision, 51: 22–45, 2015. Karsten M. Borgwardt and Hans-Peter Kriegel. Shortest-path kernels on graphs. In IEEE International Conference on Data Mining, pages 74–81, 2005. Harald Cramér and Herman Wold. Some theorems on distribution functions. Journal of the London Mathematical Society, 11(4):290–294, 1936. doi: 10.1112/jlms/s1-11.4.290. Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems, pages 2292–2300, 2013. 48
Distance-Matrix Wasserstein Statistics
Aryeh Dvoretzky, Jack Kiefer, and Jacob Wolfowitz. Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator. The Annals of Mathematical Statistics, 27(3):642–669, 1956. doi: 10.1214/aoms/1177728174. Nicolas Fournier and Arnaud Guillin. On the rate of convergence in wasserstein distance of the empirical measure. Probability Theory and Related Fields, 162:707–738, 2015. doi: 10.1007/s00440-014-0583-7. Andreas Greven, Peter Pfaffelhuber, and Anita Winter. Convergence in distribution of random metric measure spaces (Λ-coalescent measure trees). Probability Theory and Related Fields, 145(1–2):285–322, 2009. doi: 10.1007/s00440-008-0169-3. Mikhail Gromov. Metric Structures for Riemannian and Non-Riemannian Spaces. Birkhäuser, 1999. Wassily Hoeffding. Probability inequalities for sums of bounded random variables. Journal of the American Statistical Association, 58(301):13–30, 1963. doi: 10.1080/01621459. 1963.10500830. Tianyi Lin, Nhat Ho, Chenyou Fan, Marco Cuturi, and Michael I. Jordan. Projection robust wasserstein distance and riemannian optimization. In Advances in Neural Information Processing Systems, volume 33, pages 9383–9397, 2020. Pascal Massart. The tight constant in the dvoretzky–kiefer–wolfowitz inequality. The Annals of Probability, 18(3):1269–1283, 1990. doi: 10.1214/aop/1176990746. Colin McDiarmid. On the method of bounded differences. In Surveys in Combinatorics, pages 148–188. Cambridge University Press, 1989. doi: 10.1017/CBO9781107359949.008. Facundo Mémoli. Gromov–wasserstein distances and the metric approach to object matching. Foundations of Computational Mathematics, 11:417–487, 2011. Facundo Mémoli and Tom Needham. Distance distributions and inverse problems for metric measure spaces. Studies in Applied Mathematics, 149(4):943–1001, 2022. doi: 10.1111/ sapm.12526. Gabriel Peyré and Marco Cuturi. Computational optimal transport. Foundations and Trends in Machine Learning, 11(5–6):355–607, 2019. Gabriel Peyré, Marco Cuturi, and Justin Solomon. Gromov–wasserstein averaging of kernel and distance matrices. In International Conference on Machine Learning, volume 48, pages 2664–2672, 2016. I. J. Schoenberg. Metric spaces and positive definite functions. Transactions of the American Mathematical Society, 44(3):522–536, 1938. doi: 10.2307/1989894. Nino Shervashidze, S. V. N. Vishwanathan, Tobias Petri, Kurt Mehlhorn, and Karsten M. Borgwardt. Efficient graphlet kernels for large graph comparison. In International Conference on Artificial Intelligence and Statistics, pages 488–495, 2009. 49
Xu and Wu
Nino Shervashidze, Pascal Schweitzer, Erik Jan van Leeuwen, Kurt Mehlhorn, and Karsten M. Borgwardt. Weisfeiler–lehman graph kernels. Journal of Machine Learning Research, 12:2539–2561, 2011. Titouan Vayer, Laetitia Chapel, Rémi Flamary, Romain Tavenard, and Nicolas Courty. Optimal transport for structured data with application on graphs. In International Conference on Machine Learning, pages 6275–6284, 2019a. Titouan Vayer, Rémi Flamary, Romain Tavenard, Laetitia Chapel, and Nicolas Courty. Sliced gromov–wasserstein. In Advances in Neural Information Processing Systems, 2019b. Anatoly M. Vershik. Random metric spaces and universality. Russian Mathematical Surveys, 59(2):259–295, 2004. Cédric Villani. Optimal Transport: Old and New. Springer, 2009. Jonathan Weed and Francis Bach. Sharp asymptotic and finite-sample rates of convergence of empirical measures in wasserstein distance. Bernoulli, 25(4A):2620–2648, 2019. doi: 10.3150/18-BEJ1065.
50