Gromov-Wasserstein Quantization and Clustering: Structure, Rates, and Algorithms

Gromov–Wasserstein Quantization and Clustering: Structure, Rates, and Algorithms

arXiv:2608.11016v1 [math.OC] 11 Aug 2026

Florian Beier∗ and Stephan Eckstein∗

Abstract. Clustering is a fundamental class of data analysis techniques with the most important representatives being centroid-based methods like k-means. Such methods are strongly connected to quantization problems, which aim to approximate general probability measures with discrete ones. For example, k-means corresponds to quantization with respect to the Wasserstein distance. While Wasserstein quantization clusters points within a fixed space, this paper studies Gromov–Wasserstein (GW) quantization, which additionally aims at clustering the ambient geometry of the space. We show existence of solutions to the GW quantization problem and give a characterization that justifies an analogue to the k-means algorithm (Lloyd’s algorithm) to approximate them numerically. We further calculate the quantization rate for usual Euclidean geometries that are used in the GW context, and relate it to standard Wasserstein quantization rates. Finally, numerical experiments show that GW quantization opens up many modeling possibilities beyond normal clustering methods (e.g., for geodesic distances of 3D shapes or structured pruning of neural networks) and that the introduced algorithm leads to useful numerical solutions with approximation quality often in line with theoretically optimal rates.

1. Introduction. Objects in data science are increasingly large, from data sets to neural networks and high-resolution measurements. It is thus often desirable to suitably approximate these large objects with smaller ones, which are, for instance, more tractable computationally and easier to interpret. As many objects can be modeled as weighted sets of points, a common way to achieve such a reduction in size is through clustering, which replaces the points in each cluster by a single representative. This is precisely Wasserstein quantization of a finitely supported measure, for which Lloyd’s algorithm is the standard method. However, such a clustering approach, by definition, primarily cares about distances between points in a cluster and their representatives, and thus only indirectly about downstream geometric properties of the resulting object as a whole. To capture such global structure directly, we model objects as gauged measure spaces, consisting of points, weights, and a relational structure on the set of points (the gauge). The Gromov–Wasserstein (GW) distance is a type of optimal transport distance that measures the overall minimal gauge-distortion between two gauged measure spaces. This makes the GW distance well suited for comparing objects by their intrinsic relational structure. In this paper, we study the GW quantization problem for arbitrary gauged measure spaces. In contrast to the Wasserstein quantization problem, its objective is the direct and global retention of the specified gauge (i.e., the pairwise relationship between points) while reducing the object’s size. Our interest is twofold. On the practical side, moving from the Wasserstein to the GW setting affords considerable modeling freedom: an object may be encoded directly ∗

Department of Mathematics, University of Tübingen, Germany, (florian.beier@uni-tuebingen.de, stephan.eckstein@uni-tuebingen.de). Funding: FB is grateful for funding by the German Research Foundation – Project number 564105502. SE is grateful for support by the German Research Foundation through Project 553088969 as well as the Cluster of Excellence “Machine Learning — New Perspectives for Science” (EXC 2064/1 number 390727645). 1

2

F. BEIER, S. ECKSTEIN

through its pairwise relations, which need not form a metric (see Subsection 6.4 for one such example). This extends quantization to settings where ordinary clustering does not apply. For this flexibility to be practically useful, a natural goal is to lift Lloyd’s algorithm, the standard tool in the Wasserstein case, to the GW setting. On the theoretical side, it raises structural questions of independent interest, concerning the form of optimal GW quantizers and whether the optimal quantization rate matches the one in the Wasserstein setting. Main contributions. Answering the previous motivation, the main contributions of this paper are: • We introduce the GW quantization problem for general gauged measure spaces, show that it admits a solution, and characterize it in Theorem 3.4. • In Euclidean spaces of dimension d, we establish two-sided bounds relating the GW quantization value to constant multiples of the Wasserstein one (Proposition 4.1, Theorems 4.3 and 4.5), which yield the sharp n−1/d quantization rates (Corollaries 4.4 and 4.6) for approximation with n support points. • Existence of Monge maps between an input endowed with the Euclidean scalar product and its quantization is shown in Theorem 4.7. • We lift Lloyd’s algorithm to the GW setting, casting the general problem as a conditional gradient method (Algorithm 1) that, for a suitable specification, monotonically decreases the objective and has stationary cluster points (Theorem 5.3). • We validate the approach on proof-of-concept experiments: quantizing 3D shapes (Subsection 6.1) and pruning neural networks (Subsection 6.4), accelerating pairwise GW computations across many objects (Subsection 6.2), and empirically confirming the quantization rates (Subsection 6.3). Outline of the paper. The paper is structured as follows: In the remainder of this section, we discuss important related literature. In Section 2, we briefly discuss the classical Wasserstein quantization problem and its relation to clustering. In Section 3, we introduce the GW quantization problem and discuss basic properties related to existence and characterization of optimizers. Section 4 studies properties of GW quantization for Euclidean spaces and particularly the quantization rates are derived in this section. Section 5 studies the conditional gradient algorithm to (locally) solve the GW quantization problem, which is loosely inspired by Lloyd’s algorithm. Finally, Section 6 contains the numerical experiments. Supplementary statements and proofs are outsourced to the appendix. 1.1. Related Literature. Quantization of probability measures. Quantization of probability measures studies how to approximate arbitrary probability measures in Wasserstein distance with ones that have a predefined finite support size, see [30] for a monograph covering the topic. More broadly, this relates to signal processing, rate–distortion theory and vector quantization, see for instance [31] and [47]. For absolutely continuous measures on Rd with light enough tails, Zador’s theorem [77] shows that the quantization error decays at the rate n−1/d . Naturally, optimal quantization has strong relations to statistical estimation of probability measures, see for instance [16]. Closely related, the approximation of measures by point masses via discrepancies (dithering) and extensions thereof have been studied both algorithmically and through their connection to optimal transport [27, 29, 45, 66]. Finally, quantization can be restricted to

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

3

manifolds [35] or applied in general metric measure spaces [6]. We note that even in these cases, the ambient metric space is fixed and only the measure is approximated, whereas in GW quantization the ambient metric (or more generally the gauge) is approximated simultaneously. Clustering and Lloyd’s algorithm. When the measure to be approximated is already finitely supported, Wasserstein quantization coincides with k-means clustering (cf. Section 2), which is an NP-hard problem (see, e.g., [3, 39]). The most frequently employed algorithm to find a local solution is Lloyd’s algorithm [37], whose convergence is for instance studied in [51, 57]. As convergence is only local, initialization of Lloyd’s algorithm is often important [5]. Variants of the algorithm exist, like ones where cluster sizes are prescribed (corresponding to fixed weights of the approximating measure) [46], and many of the core algorithmic steps reappear in algorithms for Wasserstein barycenters [22, 73]. GW distances. Originally, GW distances were introduced by Sturm [60] and Mémoli [41] in order to analyze the space of metric measure spaces. They are distinct natural extensions of two equivalent Gromov–Hausdorff distance formulations. This paper focuses on Mémoli’s proposed GW distance which measures the overall minimal distortion required to go from one metric measure space to another. Although, in general, it amounts to numerically solving a non-convex and quadratic problem, it has been used for a wide variety of practical tasks, such as shape matching and classification [9, 41, 59], unsupervised translation [4], graph matching and node embedding [76], multi-omics [24], and particle dynamics [7]. Moreover, many theoretical insights into GW problems and variations thereof have been obtained recently [11, 23, 26, 42, 56, 68, 78]. In particular, when relaxing metrics to so-called gauges (symmetric square-integrable functions), which we also rely on in this paper, the space of gauged measure spaces endowed with the GW distance forms a complete metric space that offers a rich geometry [61]. Furthermore, as it is related to the current work, we particularly want to draw attention to the recent work on semi-discrete GW [52], where one marginal is finitely supported. Motivated by the poor tractability of solving the GW problem directly, many relaxations with faster solvers have been proposed. Examples of this are entropic regularization [28, 33, 48, 59], slicing approaches [50, 49, 69], low-rank and sample-based variants [34, 54] or iterative partitioning [75]. In the context of low-dimensional Euclidean space equipped with the squared Euclidean norm, the GW problem becomes more manageable and global solutions can be obtained in reasonable time [53]. GW quantization and related methods. GW quantization is introduced in [43] for metric measure spaces. A focus is on a type of equivalence (called duality therein) between GW quantization and clustering. Notably, both the standard GW distance from [41] and Sturm’s version [60] are treated, which lead to different results. In contrast to [43], the current work aims to establish results closer to classical Wasserstein quantization theory for general GW objectives. In particular, existence and characterization of optimal solutions, Lloyd-type iterations and general gauges are not treated in [43]. The strongest relation between the duality established in [43] and our work is given by the bounds in Section 4, where particularly some of the upper bounds could be derived from [43]. Nevertheless, while [43] focuses on bounds involving universal constants, Section 4 focuses on asymptotic rates, and importantly the main results in this section involve constants which depend on the measure which is quantized. In [67], a particular empirical GW quantization problem is studied in which all cluster centers

4

F. BEIER, S. ECKSTEIN

have the same mass. The goal of [67] is to unify both clustering and dimensionality reduction via GW quantization, which is achieved by considering GW quantization on Rp with quantized gauges being restricted to arise from points in Rd with d < p. This is in contrast to our paper, where both initial and quantized gauges are not restricted to have a euclidean structure. The quantized GW approach of [19] treats the use of any GW quantization (e.g., also suboptimal ones obtained through normal clustering) to speed up the computation of GW distances. Calculating GW barycenters [8, 10, 20], especially with fixed barycenter size [48], is a form of GW quantization (when calculating the barycenter of a single space), and some of the algorithmic steps therein relate to the ones in this paper (cf. the discussion after Theorem 3.2). Related methods further include coupling a graph to a small template, which yields graph partitioning methods [21], GW-based dictionary learning [72, 74], and the semi-relaxed GW divergence [71], which includes one free marginal as in the GW quantization problem. More broadly related approximation methods include graph coarsening [15, 18, 38, 65]. 1.2. Notation. [n] index set {1, . . . , n}, for n ∈ N. ⟨·, ·⟩, ∥ · ∥ Euclidean scalar product and norm on Rd . M+ (X) finite non-negative measures on X. P(X) probability measures on X. (2) P (X) probability measures with finite second moment (cf. Section 2). Pn (X) probability measures supported on at most n points (cf. Section 2). supp(ξ) support of ξ ∈ M+ (X), the smallest closed set of full ξ-measure. ξ|A restriction of ξ to a measurable A ⊂ X, i.e. ξ|A (B) = ξ(A ∩ B). diam(A) diameter supx,x′ ∈A dX (x, x′ ) of A ⊂ X. ξk ⇀ ξ weak convergence of finite measures. Φ# ξ push-forward of ξ by measurable Φ : X → Y , i.e. (Φ# ξ)(·) = ξ(Φ−1 (·)). PX , PY canonical projections of X × Y onto X and Y . Π(ξ, υ) couplings between ξ and υ. W2 Wasserstein distance of order 2 (cf. Section 2). GM gm-spaces up to homomorphic equivalence (cf. Section 3). (n)

Gsym

symmetric matrices G ∈ Rn×n , called gauge matrices.

GMn GW2 Πo (X, Y) Π(ξ, ∗n ) πi υi (Vi )ni=1 m(π) mξ (V ) Mξ GW qW n , qn

gm-spaces on at most n points, i.e. ([n], G, υ) with G ∈ Gsym and υ ∈ P([n]). Gromov–Wasserstein distance of order 2 (cf. Section 3). set of optimal plans between X and Y. couplings π ∈ P(X × [n]) with (PX )# π = ξ. i-th slice π(· × {i}) ∈ M+ (X) of a coupling π ∈ P(X × [n]). mass υ({i}) of υ ∈ P([n]) at i ∈ [n]. Voronoi partition of X with respect to cost functions (ci )ni=1 (cf. (3.6)). R 1 normalized mean π(X) X x dπ(x) of π ∈ M+ (X) with π(X) > 0. block mean m(ξ|V ) for measurable V ⊂ X with ξ(V ) > 0. R second-moment matrix Rd xx⊤ dξ(x) of ξ ∈ P (2) (Rd ). Wasserstein and GW quantization values ((2.1) and (3.2)).

(n)

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

QGW R

5

objective for the GW quantization problem (cf. Section 3). reduced objective inf G QGW (G, ·) of the GW quantization problem (cf. Section 5).

2. Wasserstein Quantization and Clustering. In the following, we briefly review Wasserstein quantization of probability measures and refer to [30, 6, 51] for more background. Let (X, dX ) be a Polish space, and denote by M+ (X) and P(X) the spaces of finite nonnegative and of probability measures on the Borel-σ-algebra of X induced by dX , respectively. For Polish spaces X, Y , let PX : X × Y → X and PY : X × Y → Y denote the canonical projections. For ξ ∈ P(X) and υ ∈ P(Y ), the set of couplings or transport plans of ξ and υ is Π(ξ, υ) := {π ∈ P(X × Y ) : (PX )# π = ξ, (PY )# π = υ}. The Wasserstein space (of order 2) is given by the set of measures with finite second moment Z P (2) (X) := {ξ ∈ P(X) : d2X (x0 , x) dξ(x) < ∞}, x0 ∈ X arbitrary, X

endowed with the Wasserstein distance (of order 2) Z W2 (ξ, υ) =

inf π∈Π(ξ,υ)

1

d2X (x, y) dπ(x, y)

2

,

ξ, υ ∈ P (2) (X).

X×X

In order to define the Wasserstein quantization problem, we require the set of measures on X with at most n support points, denoted by X  n n X Pn (X) := υi δyi : yi ∈ X, υi ≥ 0, i = 1, . . . , n, υi = 1 . i=1

i=1

Let n ∈ N and ξ ∈ P (2) (X), the Wasserstein quantization problem reads (2.1)

qW n (ξ) :=

inf υ∈Pn (X)

W2 (ξ, υ).

While more general distances can be employed, in this paper we focus on just W2 (and later the GW distance of order 2). For convex X ⊂ Rd , existence of quantizations, i.e. solutions υ ∗ 1 PN to (2.1), is guaranteed by [30, Thm. 4.12]. If further ξ = N i=1 δxi ∈ P(X) is an empirical measure, the quantization problem (2.1) is precisely Euclidean k-means (minimum-sum-ofsquares clustering) with k = n, see [30], which is NP-hard [3]. P The interest of solving (2.1) is twofold. Firstly, a solution υ = ni=1 υi δyi ∈ Pn (X) provides a quantization (i.e. a smaller approximate representation) of ξ according to the Wasserstein distance W2 . Secondly, any π ∈ Π(ξ, υ), that solves the inner problem provides a clustering of ξ according to dX . More precisely, the measures πi := π(· × {yi }) describe the part of the measure ξ which belongs to the i-th cluster with center yi . With this in mind, it becomes clear that we can rewrite n Z X W 2 (2.2) qn (ξ) = inf d2X (x, yi ) dπi (x). n (yi )i=1 ⊂X Pn i=1 (πi )n i=1 ⊂M+ (X), i=1 πi =ξ

X

6

F. BEIER, S. ECKSTEIN

Minimizing out the (πi )ni=1 by assigning each point to a nearest center, (2.2) can be reduced to the n-point covering form Z W 2 (2.3) qn (ξ) = inf d2X (x, Y ) dξ(x), dX (x, Y ) := min dX (x, y). Y ⊂X X |Y |≤n

y∈Y

Lloyd’s algorithm [37] is based on an alternating optimization of (yi )ni=1 and (πi )ni=1 . Optimizing over one parameter and keeping the other one fixed allows for closed form solutions. The first insight in this regard is that, given (yi )ni=1 , πi is concentrated on {x ∈ X : dX (x, yi ) = minj=1,...,n dX (x, yj )}, which gives rise to the importance of so-called Voronoi partitions of the space X. Conversely, the second insight is that given πi , the yi are simply given as suitable means of the parts of ξ corresponding to πi . For X ⊂ Rd convex and the Euclidean distance dX (·, ·) = ∥ · − · ∥, this is just the standard arithmetic mean, i.e. (2.2) simplifies further to (2.4)

2 qW n (ξ) =

inf

n Z X

V1 ,...,V Snn ⊂X disjoint i=1 i=1 Vi =X

2

∥x − mξ (Vi )∥ dξ(x),

mξ (Vi ) :=

Vi

1 ξ(Vi )

Z x dξ(x) Vi

with the convention that mξ (Vi ) ∈ X is set arbitrarily for ξ(Vi ) = 0. We derive analogues to these properties for GW quantization in Section 3, which are the basis for the algorithmic ideas in Section 5. Regarding approximation quality, Zador’s theorem characterizes the asymptotic rate of approximation, i.e. behavior of the value qW n (ξ) for n → ∞. More precisely, consider R the2+δ d ξ ∈ P(R ) with ∥x∥ dξ(x) < ∞ for some δ > 0 and assume that ξ contains an absolutely continuous part with non-zero mass. Then, Zador’s theorem [30, Thm. 6.2] yields limn→∞ n1/d qW n (ξ) ∈ (0, ∞). We exploit this in Section 4 to show similar rates of convergence for GW quantization problems in common Euclidean settings as well. 3. Gromov–Wasserstein Quantizations: Existence and Characterization . The GW distance [41, 61] is a type of optimal transport problem which aligns and compares based on intrinsic (dis-)similarities as opposed to a given ambient metric. After giving a brief summary and fixing the notation, we formulate and analyze the GW quantization problem below. A gauged measure space (gm-space) is a triple (X, g, ξ), where X is a Polish space, ξ ∈ P(X) is a probability measure and g ∈ L2sym (ξ ⊗ ξ) is an associated gauge, where L2sym (ξ ⊗ ξ) := {g ∈ L2 (ξ ⊗ ξ) : g is symmetric}. We denote the set of gm-spaces by GM. For two gm-spaces X = (X, g, ξ) ∈ GM, Y = (Y, h, υ) ∈ GM, the GW distance (of order 2) is given as Z Z 1 2 (3.1) GW2 (X, Y) := inf |g(x, x′ ) − h(y, y ′ )|2 dπ(x, y) dπ(x′ , y ′ ) . π∈Π(ξ,υ)

(X×Y )2

The infimum in (3.1) is always attained and the set of solutions is denoted as Πo (X, Y). The GW distance is only a pseudometric on gm-spaces. More precisely, GW2 (X, Y) may vanish for distinct X = (X, g, ξ) and Y = (Y, h, υ). The gm-spaces X, Y with GW2 (X, Y) = 0

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

7

are called homomorphic which admits an exact characterization, for which we refer to [61]. For our purposes only the following sufficient condition is needed. If there exists a Borelmeasurable map Φ : X → Y with υ = Φ# ξ and h(Φ(x), Φ(x′ )) = g(x, x′ ) for ξ ⊗ ξ-almost every (x, x′ ), then X and Y are homomorphic. We call such a map Φ a homomorphism from X to Y. By identifying homomorphic gm-spaces, GW2 becomes a complete metric on the resulting quotient of GM. Throughout, we write X ∈ GM for a gm-space as well as for its (homomorphism) equivalence class. Let n ∈ N and set [n] := {1, . . . , n}. Following [43] for metric measure spaces, we want to consider the natural generalization of the Wasserstein quantization problem (2.1) to the GW case for which we consider gm-spaces Y = (Y, h, υ) ∈ GM whose space Y has n ∈ N elements. As renaming the elements constitutes a homomorphism, any such gm-space Y always admits a representative whose space is [n]. Moreover, we identify the gauges h for such spaces as gauge matrices (n) := {G ∈ Rn×n : G symmetric}, G ∈ Gsym

via h(i, i′ ) := Gi,i′ , i, i′ ∈ [n]. The space of gm-spaces that are supported on at most n points GMn ⊂ GM is given by  (n) GMn := ([n], G, υ) : G ∈ Gsym , υ ∈ P([n]) . Let n ∈ N and X ∈ GM be an arbitrary gm-space. We consider the GW quantization problem

(3.2)

qGW n (X) :=

inf GW2 (X, Y).

Y∈GMn

Since GW is invariant under homomorphisms, so is qGW n . Before we show existence of solutions and characterize them, we make the following observation in terms of the support size of elements in GMn . Proposition 3.1. Let Y ∈ GMn . There exists a representative Ỹ = ([n], G̃, υ̃) of Y with fully supported measure υ̃ ∈ P([n]). The proof, presented in Appendix B.1, repeatedly duplicates a positive-mass point of Y into two copies that inherit its gauge values and share its mass. Each such split leaves the gm-space unchanged up to homomorphism while enlarging the support, until all n points carry positive mass. Proposition 3.1 allows us to always consider fully supported representatives of any elements in GMn which we will do throughout. Similar to Section 2, we encounter transport plans of the (semi-discrete) form π ∈ P(X × [n]), where X is some Polish space. In this setting, we use the notation πi := π(· × {i}) ∈ M+ (X), as well as υi := υ({i}) for any υ ∈ P([n]), i ∈ [n]. With this notation, the marginal of π on [n] can be expressed as (P[n] )# π = π(X × ·) =

n X i=1

π(X × {i}) · δi =

n X i=1

πi (X) · δi .

8

F. BEIER, S. ECKSTEIN

Furthermore, the marginal on X can be expressed as (PX )# π = π(· × [n]) =

n X

πi .

i=1

For ξ ∈ P(X), we set Π(ξ, ∗n ) := {π ∈ P(X × [n]) : (PX )# π = ξ}. Since the discrete space is fixed as [n], the measure positions of (P[n] )# π, π ∈ Π(ξ, ∗n ), do not vary. Instead their geometry varies via the associated Pn gauge matrix. This allows us to uniquely identify n π ∈ Π(ξ, ∗n ) by (πi )i=1 ⊂ M+ (X), i=1 πi = ξ. This stands in contrast to the Wasserstein setting, where the measure positions vary directly. Throughout, we will frequently make use of this one-to-one correspondence in the GW setting, that is, we frequently consider π = (πi )ni=1 . We turn our attention to showing existence of GW quantizations. We follow a similar alternating approach as the one that Lloyd’s algorithm for the Wasserstein case is based on. That is, we analyze the minimization with respect to the gauge matrix and the inner GW transport in the quantization problem (3.2) disjointly at first. For X = (X, g, ξ), we set n ZZ X n QGW (G, (πi )i=1 ) := (g(x, x′ ) − Gi,i′ )2 dπi (x) dπi′ (x′ ), i,i′ =1

X2

so that the GW quantization problem (3.2) can be written as (3.3)

2 qGW n (X) =

QGW (G, (πi )ni=1 ).

inf (n)

G∈GsymP n (πi )n i=1 ⊂M+ (X), i=1 πi =ξ

In the Wasserstein case, the minimization with respect to the measure positions admits closed form. Analogously for the GW case, we obtain a similar result for the minimization with respect to the gauge matrix. Theorem 3.2. Let X = (X, g, ξ) ∈ GM be arbitrary and fix π ∈ P(X ×[n]) with (PX )# π = ξ and πi (X) > 0 for all i ∈ [n]. Then inf G∈G (n) QGW (G, (πi )ni=1 ) is uniquely solved by sym

(3.4)

(n) Ĝ ∈ Gsym ,

1 Ĝi,i′ := πi (X)πi′ (X)

ZZ

g(x, x′ ) dπi (x) dπi′ (x′ ),

i, i′ ∈ [n].

X2

(n)

Set Y = ([n], G, υ) for G ∈ Gsym arbitrary and Ŷ := ([n], Ĝ, υ), where υ = (P[n] )# π. If additionally π ∈ Πo (X, Y), then the following estimate holds (3.5)

GW22 (X, Y) ≥ GW22 (X, Ŷ) + ∥G − Ĝ∥2L2 (υ⊗υ) .

The proof, presented in Appendix B.2, hinges analogously to the Wasserstein setting on the identification of the L2 mean. We note that the optimal gauge Ĝ coincides with the barycenter-matrix update of the GW barycenter algorithm in [48] for a single input space and that the same formula also occurs in [18] in a framework related to graph coarsening. For cost functions (ci )ni=1 on X, the associated Voronoi partition (Vi )ni=1 of X is defined as  (3.6) Vi := T −1 ({i}), T (x) := min i ∈ [n] : ci (x) = minj∈[n] cj (x) .

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

9

In contrast to the Wasserstein case, the minimization of QGW (G, ·) over Π(ξ, ∗n ) for a fixed (n) gauge matrix G ∈ Gsym does not admit closed form. This is due to the fact that the problem is quadratic and non-convex. Nevertheless, we can show the following characterization. (n)

Theorem 3.3. Let X = (X, g, ξ) ∈ GM be arbitrary and G ∈ G Psym be fixed. Then QGW (G, ·) attains its minimum over Π(ξ, ∗n ). Let (π̂i )ni=1 ⊂ M+ (X) with ni=1 π̂i = ξ be any minimizer. Then, for every i ∈ [n], the measure π̂i is concentrated on (π̂)

(π̂)

Ai := {x ∈ X : ci (x) = min cj (x)}, j

(π̂)

where ci (x) :=

n Z X

(g(x, x′ ) − Gi,i′ )2 dπ̂i′ (x′ ).

i′ =1 X

If additionally ξ(Ai ∩ Aj ) = 0 for all i ̸= j, then π̂i = ξ|Vi , i ∈ [n], where (Vi )ni=1 is the Voronoi (π̂) partition of X with respect to (ci )ni=1 . The proof is given in Appendix B.3. Existence follows from standard arguments. For the characterization, the key step is to show that a minimizer π̂ of the quadratic problem also minimizes its linearization at π̂, whose minimizers satisfy the given properties. We remark that the condition ξ(Ai ∩ Aj ) = 0 cannot be dropped. If Gi,i′ = c for all ′ i, i ∈ [n], then QGW (G, (πi )ni=1 ) =

n ZZ X i,i′ =1

′

ZZ

2

(g(x, x′ ) − c)2 dξ dξ

(g(x, x ) − c) dπi dπi′ =

X2

X2

is independent of (πi ). Thus, every feasible coupling is a minimizer, and correspondingly, R (π̂) ci ≡ X (g(·, x′ ) − c)2 dξ is independent of i so that Ai = X, i ∈ [n], and ξ(Ai ∩ Aj ) = 1, i ̸= j. In particular, minimizers need not be Voronoi partitions. On the basis of the previous analysis, we are now able to establish existence for the GW quantization problem (3.2) and give a characterization of solutions. Theorem 3.4. Let X ∈ GM. The quantization problem (3.2) admits a solution. Moreover, let Ŷ = ([n], Ĝ, υ̂) ∈ GMn be any solution of (3.2), with fully supported υ̂ ∈ P([n]), see Proposition 3.1, Z Z and let π̂ ∈ Πo (X, Ŷ). Then the following hold: 1 (i) Ĝi,i′ = g(x, x′ ) dπ̂i (x) dπ̂i′ (x′ ), i, i′ ∈ [n]. υ̂i υ̂i′ 2 X (ii) For every i ∈ [n], the measure π̂i is concentrated on n o Ai := x ∈ X : ĉi (x) = min ĉj (x) , j∈[n]

2 2 (iii) qGW n (X) = GW2 (X, Ŷ) =

ZZ X2

where

ĉi (x) :=

n Z X

g(x, x′ ) − Ĝi,i′

2

dπ̂i′ (x′ ).

i′ =1 X

g 2 (x, x′ ) dξ(x) dξ(x′ ) −

n X

Ĝ2i,i′ υ̂i υ̂i′ .

i,i′ =1

If additionally ξ(Ai ∩ Aj ) = 0 for all i ̸= j, then π̂ = (Id, T )# ξ, where T (x)P = i if x ∈ Vi , for n n the Voronoi partition (Vi )i=1 of X with respect to (ĉi )i=1 . In this case, υ̂ = ni=1 ξ(Vi ) · δi and ZZ 1 Ĝi,i′ = g(x, x′ ) dξ(x) dξ(x′ ), i, i′ ∈ [n]. ξ(Vi )ξ(Vi′ ) Vi ×Vi′

10

F. BEIER, S. ECKSTEIN

The strategy for the proof, presented in Appendix B.4, is as follows. General existence of GW quantizers is shown by taking a minimizing sequence and applying the pointwise bound from Theorem 3.2. Then convergence of the gm-spaces reduces to convergence of the measures which subsequently holds due to weak compactness. The characterization is obtained by using the connection of the quantization problem with QGW and the established properties of pointwise minimizers of the latter, i.e. Theorems 3.2 and 3.3. The characterization in Theorem 3.4 induces a self-consistency condition. More precisely, each π̂i is concentrated on the cells Ai determined by (ĉi )ni=1 , while the costs ĉi in turn depend on π̂. Optimal couplings are therefore fixed points of the assignment map sending a coupling to the one whose slices are concentrated on Ai , i ∈ [n]. In general, not every such fixed point is a global minimizer, and identifying the spurious ones is to our knowledge an open problem. A fixed-point characterization of the same flavor underlies the computation of the tangential GW barycenters in [8]. 4. Gromov–Wasserstein Quantization for Euclidean Spaces. We turn our attention to Euclidean settings, that is, we consider the quantization problem (3.2) for X = (X, g, ξ) ∈ GM, where X ⊂ Rd , ξ ∈ P(Rd ) and g is either the Euclidean distance, the squared S Euclidean distance or the Euclidean scalar product. The set of finite gm-spaces, i.e. ∞ n=1 GMn , is GW dense in GM [61, Thm. 5.28]. So for the quantization value qn (X) defined in (3.2), it holds qGW n (X) → 0, as n → ∞, for every X ∈ GM. Below we quantify this convergence, i.e. we show that under weak assumptions on X ⊂ Rd and ξ ∈ P(Rd ), qGW n (X) behaves asymptotically like 1 n− d , for the aforementioned Euclidean gauges. In each of these cases, we do this by showing the existence of constants c, C > 0 for which (4.1)

GW W c · qW n (ξ) ≤ qn (X) ≤ C · qn (ξ).

The constant for the lower bound c = c(ξ) hereby always depends on ξ (and additionally on the diameter of X for the Euclidean gauge), and the constant C is universal for the Euclidean gauge, depending on X for the squared Euclidean gauge, and depending on ξ for the scalar product gauge. Here qW n (ξ) is taken with respect to the Euclidean metric in every case. The asymptotic rate of qGW n (X) is then established by Zador’s theorem ([30, Thm. 6.2]) −1 which ensures that generally qW n (ξ) behaves asymptotically like n d . For the special case of the scalar product, we also establish the existence of Monge maps from X to solutions of the quantization problem (3.2) and give a one-to-one correspondence to the Wasserstein quantization problem in one dimension. All proofs of this section are given in Appendix C. 4.1. Quantization Rate for the Euclidean Distance and Its Square. In the following we establish the two-sided bounds (4.1) for g ∈ {∥ · − · ∥, ∥ · − · ∥2 }. The upper bounds rely on the following more general result on Lipschitz transformations of metrics. Proposition 4.1. Let X = (X, g, ξ) ∈ GM be a gm-space with g = φ ◦ dX , where dX is a metric on X and φ : [0, ∞) → R. Let qW n (ξ) be the Wasserstein quantization value defined in (2.1) with respect to the metric dX . If φ is L-Lipschitz continuous on the range of dX , then (4.2)

W qGW n (X) ≤ 2L qn (ξ).

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

11

The proof relies on applying the Lipschitz bound together with the triangle inequality to bound the double GW-like integral from above by a single Wasserstein-like integral. It is presented in Appendix C.1. Remark 4.2. Consider X = (X, ∥ · − · ∥p , ξ) ∈ GM for some p ≥ 1. If X ⊂ Rd is bounded, W the previous result implies an upper bound of qGW n (X) via qn (ξ) with respect to the standard Euclidean metric. More precisely, on [0, diam(X)], the function φ(a) = ap satisfies |φ′ (a)| = p ap−1 ≤ p diam(X) p−1 , and is therefore p diam(X) p−1 -Lipschitz continuous. In this case, we obtain p−1 W qGW qn (ξ). n (X) ≤ 2p diam(X) W GW W In particular, qGW n (X) ≤ 2 qn (ξ) for p = 1 and qn (X) ≤ 4 diam(X) qn (ξ) for p = 2.

In the following, λmin (M ) and λmax (M ) denote the smallest, respectively largest, eigenvalue of a positive semi-definite matrix M . We turn our attention to the derivation of the lower bounds. Theorem 4.3. Let X = (X, g, ξ) ∈ GM with X ⊂ Rd closed convex. The mean and covariR R ance matrix of ξ ∈ P(X) are set as m̄ξ := X x dξ(x), and Σξ := X (x − m̄ξ )(x − m̄ξ )⊤ dξ(x), respectively. Let qW n (ξ) be defined as in (2.1) with respect to the standard Euclidean metric. R 4 (i) If X ∥x∥ dξ(x) < ∞ and g = ∥· − ·∥2 , then q λmin (Σξ ) · qW qGW (X) ≥ 2 n n (ξ). (ii) If X is bounded and g = ∥ · − · ∥, then qGW n (X) ≥

p λmin (Σξ ) W · qn (ξ). diam(X)

The proof, presented in Appendix C.2, relies on lower bounding the GW value as the average within-cell variance of the gauge g(·, y), y ∈ X. The latter can in turn be related to the Wasserstein quantization value. Corollary 4.4. Let X = (X, g, ξ) ∈ GM with X ⊂ Rd bounded and ξ ∈ P(X) be absolutely −1/d , that is continuous. For g ∈ {∥ · − · ∥, ∥ · − · ∥2 }, it holds qGW n (X) ≍ n lim sup n1/d qGW n (X) < ∞ n→∞

and

lim inf n1/d qGW n (X) > 0. n→∞

Proof. Without loss of generality, we assume that X ⊂ Rd is additionally closed and convex. Otherwise, we can exchange X with the representative whose base space is the closed convex hull of X and g, ξ unchanged. Firstly, note that the convergence of qW n (ξ) with respect − d1 d to the standard Euclidean metric on R is n due to Zador’s theorem ([30, Thm. 6.2]). So it suffices in both cases to bound qGW (X) from below and above by qW n n (ξ) up to a positive multiplicative constant. The upper bound follows directly by Proposition 4.1, see Remark 4.2. The lower bounds follow by Theorem 4.3. It remains to establish that the smallest eigenvalue

12

F. BEIER, S. ECKSTEIN

R R of Σξ = X (x − m̄ξ )(x − m̄ξ )⊤ dξ(x), with m̄ξ = X x dξ(x), is non-zero. For any v ∈ Rd with v ̸= 0, it holds Z (v ⊤ (x − m̄ξ ))2 dξ(x) > 0.

v ⊤ Σξ v =

X

The last inequality holds as {x ∈ X : v ⊤ (x − m̄ξ ) = 0} is a Lebesgue null-set. Therefore, λmin (Σξ ) > 0 as desired. 4.2. The Scalar Product: Quantization Rate and Monge Maps. We turn our attention to the scalar product case, i.e. throughout we consider X = (X, ⟨·, ·⟩, ξ) ∈ GM with X ⊂ Rd . In this setting, the quantization problem admits an accessible alternative formulation. It allows us to derive two-sided bounds of the form (4.1) and establish the existence of Monge maps from X to a quantizer. Theorem 4.5. Let X ⊂ Rd be closed convexR and X = (X, ⟨·, ·⟩, ξ) ∈ GM. Consider the symmetric positive semi-definite matrix Mξ = X xx⊤ dξ(x) ∈ Rd×d . The GW quantization value admits the bounds q q W GW λmin (Mξ ) · qn (ξ) ≤ qn (X) ≤ 2λmax (Mξ ) · qW n (ξ), where qW n (ξ) is defined in (2.1) with respect to the standard Euclidean metric. The proof of this result is provided in Appendix C.3. Therein, we also establish that ξ ∈ P (2) (X), so that Mξ and qW n (ξ) are well-defined. The proof relies on auxiliary results 2 as an infimum of the within-cell variance which is weighted by a matrix (X) that express qGW n 2 contained between Mξ and 2Mξ . This bounds qGW n (X) between λmin (Mξ ) and 2λmax (Mξ ) 2 times the unweighted within-cell variance, whose infimum over couplings is qW n (ξ) . Corollary 4.6 (Quantization Let X = (X, ⟨·, ·⟩, ξ) ∈ GM with X ⊂ Rd and ξ ∈ P(X) R rate). −1/d , that is 2+δ be absolutely continuous. If X ∥x∥ dξ < ∞ for some δ > 0, then qGW n (X) ≍ n lim sup n1/d qGW n (X) < ∞ n→∞

and

lim inf n1/d qGW n (X) > 0. n→∞

Proof. Due to the bounds from Theorem 4.5, the proof is analogous to Corollary 4.4. The lower constant in (4.1) is necessarily measure-dependent. In Appendix E we exhibit a family of gauged measure spaces X(n) with measures ξ (n) , valid for all three Euclidean gauges, (n) )/qW (ξ (n) ) → 0, ruling out any uniform lower constant c > 0. with qGW n (X n The existence of Monge maps—optimal couplings concentrated on the graph of a map—is a recurring theme in the realm of optimal transport. In the GW setting, for the inner-product cost, [26] establish the existence of an optimal Monge map under absolute continuity of the source measure, whereas for the squared-distance cost they obtain, in general, only an optimal 2-map. The latter is a plan that couples source masses to at most 2 targets. In the semidiscrete, inner-product setting, [52] likewise establish the existence of Monge maps, but under absolute continuity of ξ together with a non-degeneracy condition on the discrete marginal (guaranteed, e.g., when its support points have pairwise distinct norms). In the following, we present a result ensuring the existence of Monge maps for the quantization problem (3.2)

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

13

in the scalar-product setting that requires neither absolute continuity of ξ nor any such nondegeneracy condition, since here the discrete object is itself subject to optimization. For the following, recall Z 1 mξ (V ) = x dξ(x), V ⊂ X measurable, ξ(V ) V where mξ (V ) ∈ X is set arbitrarily for ξ(V ) = 0. Theorem 4.7 (Monge Maps for Scalar Product Spaces). Let X = (X, ⟨·, ·⟩, ξ) ∈ GM with X ⊂ Rd . Then there exists a measurable map T : X → [n] such that, for the disjoint partition Vi := T −1 ({i}), i ∈ [n], the gm-space Ŷ := ([n], Ĝ, υ̂) ∈ GMn ,

Ĝ

i,i′

:= ⟨mξ (Vi ), mξ (Vi′ )⟩,

υ̂ :=

n X

ξ(Vi ) · δi ,

i=1

is a solution of the GW quantization problem (3.2) with respect to X, and (Id, T )# ξ ∈ Πo (X, Ŷ). The proof, presented in Appendix C.4, establishes the result by exploiting the reformulation of qGW n (X) that we established in Appendix C.3 for the proof of Theorem 4.5. The latter characterizes qGW as the minimization of a concave function over the set Π(ξ, ∗n ). By n the Bauer maximum principle, an extreme point of Π(ξ, ∗n ) is a solution. Showing that the extreme points are map-induced plans is a standard argument. The previous theorem shows that in the scalar-product setting X = (X, ⟨·, ·⟩, ξ) with X ⊂ Rd , an optimal quantizer is induced by a partition (Vi )ni=1 of X: the optimal coupling is the hard assignment πi = ξ|Vi and the optimal gauge is the block-mean gauge Ĝi,i′ = ⟨mξ (Vi ), mξ (Vi′ )⟩, i, i′ ∈ [n]. The quantization problem therefore reduces to the partition problem (4.3)

2 qGW n (X) =

inf

n Z X

V1 ,...,V Snn ⊂X disjoint i,i′ =1 Vi ×Vi′ i=1 Vi =X

2 ⟨x, x′ ⟩ − ⟨mξ (Vi ), mξ (Vi′ )⟩ dξ(x) dξ(x′ ).

This is the natural GW analogue of the Wasserstein partition formulation (2.4): the Wasserstein objective penalizes the squared deviation of each point x from its cell barycenter mξ (Vi ), whereas the GW objective penalizes the squared deviation of each pairwise scalar product value ⟨x, x′ ⟩ from the scalar product value of the block means ⟨mξ (Vi ), mξ (Vi′ )⟩ over Vi × Vi′ . We conclude the section by noting that, in one dimension, the Wasserstein and GW quantization problems reduce to the same problem. Theorem 4.8. Let X = (X, ⟨·, ·⟩, ξ) ∈ GM with X ⊆ span(u), u ∈ Rd , ∥u∥ = 1, and set R T : X → R, T (x) = u⊤ x as well as M2 (T# ξ) := R t2 d(T# ξ)(t). Then
2 W 2 W 2 qGW n (X) = qn (T# ξ) 2M2 (T# ξ) − qn (T# ξ) . The proof, presented in Appendix C.5, uses that X is homomorphic to the one-dimensional gm-space (R, (s, t) 7→ st, T# ξ) and a direct reformulation of the appearing integrals in one dimension.

14

F. BEIER, S. ECKSTEIN

Algorithm 1 Conditional Gradient for GW Quantization 1: Input: gm-space X := (X, g, ξ) number of support points n ∈ N 2: Initialize π ∈ Π(ξ, ∗n ) with πi (X) > 0 and set υi ← πi (X) 3: while not converged do RR 4: Ĝ(π)i,i′ ← υi1υ ′ X 2 g(x, x′ ) dπi (x) dπi′ (x′ ), ▷ indices i, i′ range over supp(υ) throughout i
R P 2 (π) 5: ĉi ← ni′ =1 X g(·, x′ ) − Ĝ(π)i,i′ dπi′ (x′ ), (π)

πi ← (1 − t) πi + t ξ|Vi for a Voronoi partition (Vi )ni=1 of X with respect to (ĉi )ni=1 and a step size t ∈ [0, 1], 7: υi ← πi (X) 8: end while 9: Output: Quantizer Y = ([n], Ĝ(π), υ), coupling π ∈ Π(ξ, υ) ▷ Ĝ(π) = 0 outside supp(υ)2

6:

5. Algorithm and Convergence. We turn to the practical approximation of solutions to the GW quantization problem (3.2) for which we propose Algorithm 1. While the initial idea for the algorithm arose as an analogue to Lloyd’s algorithm, we formally derive the algorithm as a conditional gradient method of a reduced objective R of QGW defined below. The algorithm is stated for generic convergence criteria. For our experiments, we observe the Voronoi partitions becoming stationary fairly quickly, which we use as the convergence criterion. Regarding the choice of the step size, see the discussion after Proposition 5.1. Consider an input gm-space X := (X, g, ξ) and fix the support size of the quantization n ∈ N. We define the gauge-reduced objective (5.1)

R(π) :=

inf (n)

π ∈ Π(ξ, ∗n ).

QGW (G, π),

G∈Gsym

For every π ∈ Π(ξ, ∗n ) with πi (X) > 0, i ∈ [n], the minimization in (5.1) is uniquely attained at ZZ 1 Ĝ(π)i,i′ := g dπi dπi′ , i, i′ ∈ [n], πi (X) πi′ (X) X 2
2 see Theorem 3.2. Thus, R(π) = QGW Ĝ(π), (πi )ni=1 and qGW n (X) = inf π∈Π(ξ,∗n ) R(π). To analyze the gradient of R, set n Z X
2 (π) ĉi (x) := g(x, x′ ) − Ĝ(π)i,i′ dπi′ (x′ ), i ∈ [n]. i′ =1 X

(π)

Since πi′ ≤ ξ, i′ ∈ [n], and g ∈ L2 (ξ ⊗ ξ), we have ĉi ∈ L1 (ξ). The directional derivative of the reduced objective R can be derived similarly to the one of the standard GW objective [48, 68]. Proposition 5.1. Let X = (X, g, ξ) ∈ GM, n ∈ N, and let π ∈ Π(ξ, ∗n ) with πi (X) > 0 for all i ∈ [n]. Then, for every η ∈ Π(ξ, ∗n ), the one-sided directional derivative of R at π exists and equals
n Z X R (1 − t)π + tη − R(π) (π) lim =2 ĉi d(ηi − πi ). t↓0 t X i=1

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

15

The proof is given in Appendix D.1. It relies on the fact that, along the segment πt = (1 − t)π + tη, the block-mean gauge Ĝ(πt ) is a ratio of polynomials of t with non-vanishing denominator, and thus differentiable. Differentiating R at t = 0 and simplifying yields the stated linear form. We minimize R over the convex set Π(ξ, ∗n ) by conditional gradient descent: at a current π we move toward a minimizer of the linearization of R, that is, of the directional derivative Pn R (π) of Proposition 5.1. The linear functional η 7→ i=1 X ĉi dηi is minimized over Π(ξ, ∗n ) (π)

(π)

exactly by those η, such that ηi is concentrated on {x ∈ X : ĉi (x) = minj ĉj (x)} in particular by the Voronoi assignment (ξ|Vi )ni=1 (cf. Lemma A.1), i ∈ [n]. Hence Algorithm 1 is the conditional-gradient method for minπ∈Π(ξ,∗n ) R(π): line 4 evaluates the gradient through Ĝ(π), lines 5–6 determine a minimizer ξ|Vi of the linearization, and the convex combination in line 6 is the conditional-gradient step. Taking the unit step size (t = 1) in line 6 ensures the hard assignment πi = ξ|Vi which may be of particular interest when seeking to cluster X. In this case, it may happen that clusters attain zero mass, i.e. υi = 0. This contrasts the case t ∈ [0, 1) where clusters masses remain positive throughout. We also note, that the unit step t = 1 need not decrease R = QGW (Ĝ(·), ·) in general. More precisely, by Proposition 5.1 the conditional-gradient direction is one of descent, so R decreases for small t > 0, but without concavity the full step may increase it. Remark 5.2 (Line Search). Let (Vi )ni=1 be the Voronoi partition of X with respect to (π) (ĉi )ni=1 . A monotone decrease of R is guaranteed by the line search
t ∈ argmin QGW Ĝ(π), π + s(η − π) ,

ηi := ξ|Vi ,

i ∈ [n].

s∈[0,1]

This is the minimization of the quadratic functional s 7→ R(π) − s b + s2 a, where b := −2

n Z X i=1

X

(π) ĉi d(ηi − πi ),

a :=

n ZZ X i,i′ =1

g − Ĝ(π)i,i′

2

d(ηi − πi ) d(ηi′ − πi′ ).

X2

Note that b ≥ 0 since η minimizes the linearization.
Hence t = 1 if a ≤ 0 and t = min{1, b/(2a)} if a > 0, and in either case R (1 − t)π + tη ≤ R(π). We can use the line search as described in Remark 5.2 to produce a sequence (π (k) )k∈N which monotonically decreases R, that is R(π (k+1) ) ≤ R(π (k) ). Since R is non-negative, this ensures convergence of (R(π (k) ))k∈N . At the same time, since Π(ξ, ∗n ) is weakly compact, the sequence (π (k) )k∈N admits cluster points. Under certain assumptions, these cluster points are also stationary for R as the following result shows. Theorem 5.3. Let X = (X, g, ξ) ∈ GM with bounded gauge g. Let n ∈ N, and let (π (k) )k ⊂ Π(ξ, ∗n ) be generated by Algorithm 1 with line search as described in Remark 5.2. Assume the

16

F. BEIER, S. ECKSTEIN (k)

iterates satisfy πi (X) > 0 for all i ∈ [n], k ∈ N. Then every cluster point π ∗ of (π (k) )k∈N with πi∗ (X) > 0 for all i ∈ [n] is stationary, that is, no feasible direction is one of descent: 2

n Z X i=1

X

(π ∗ )

ĉi

d(ηi − πi∗ ) ≥ 0

for all η ∈ Π(ξ, ∗n ).

The proof is given in Appendix D.2 and follows the standard conditional-gradient (Frank– Wolfe) argument. More precisely, the line search enforces monotone decrease of R, so by Proposition 5.1 no feasible direction at a cluster point is one of descent. Note that the (k) positivity assumption πi (X) > 0, i ∈ [n], k ∈ N is satisfied whenever the step sizes are strictly less than 1 when starting from an initial π ∈ Π(ξ, ∗n ) with πi (X) > 0, for all i ∈ [n]. We remarked before that for the unit step size, monotone decrease of R is not ensured. The following result shows that the scalar product case is an exception, where concavity of R renders the unit step monotone. Corollary 5.4. Let X ⊂ Rd and g = ⟨·, ·⟩. Then R is concave on Π(ξ, ∗n ). Let (π (k) )k∈N be (k) produced by Algorithm 1 with unit step size t = 1, and assume πi (X) > 0 for all i ∈ [n], k ∈ N. Then R(π (k+1) ) ≤ R(π (k) ), k ∈ N. If moreover X is finite, the iteration reaches, in finitely many steps, a coupling π ⋆ such that, (π ⋆ ) (π ⋆ ) for each i ∈ [n], πi⋆ is concentrated on {x ∈ X : ĉi (x) = minj ĉj (x)}. The proof is given in Appendix D.3. Concavity bounds the unit-step increment by the directional derivative of Proposition 5.1, which is non-positive since the update minimizes the linearization (Lemma A.1), giving monotone decrease. For finite X the iterates are polytope vertices across which R strictly decreases, forcing termination. 6. Numerical Examples. In the following, we provide various numerical experiments. These experiments serve primarily three goals: First, showcase that the proposed GW quantization is numerically feasible and leads to structurally different results compared to normal Wasserstein clustering (cf. particularly Subsection 6.1). Second, show that the added modeling freedom provided by gm-spaces is useful for certain applications and that the quantization problem can also be solved effectively in such settings (cf. Subsection 6.4). And third, show that in certain simple cases, the proposed algorithm leads to approximation rates in line with the theory for optimal quantization (cf. Subsection 6.3) and that the approximation quality leads to desirable benefits in downstream tasks (cf. Subsection 6.2). We emphasize that all experiments should be considered as proof-of-concept, and not (yet) improving on state-of-the-art performance for important practical tasks. Nevertheless, we believe these experiments nicely showcase the practical possibilities that GW quantization can provide for data science applications. All experiments are implemented in Python.1 In addition to the self-implemented Algorithm 1, we use the fast simplification package for mesh decimation, as well as POT [28] for 1

The source code is publicly available at https://github.com/florian-beier/gw-quantization.

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

17

some GW distance and coupling computations. Accordingly, the reported values for the GWproblem are numerical approximations of the true values arising from the respective algorithms used in the POT package. 6.1. Quantization of a 3D Shape. In this example, we quantize a 3D surface and visualize the resulting quantizers at decreasing resolutions. The surface is a 3D registration of a human subject taken from the training part of the FAUST dataset [13]. It consists of 6890 vertices and 13 776 faces. We regard the surface as a graph whose nodes are the mesh vertices and whose edges are the edges of the triangulation. Each edge is weighted by the Euclidean distance between the vertices it connects. We convert the 3D shape to a gm-space X = (X, g, ξ) as follows. The space X ⊂ R3 is set as the vertex set, the gauge g consists of the pairwise Dijkstra distances on the weighted graph, and ξ is the normalized area-weighted vertex measure, i.e. the mass of a vertex is proportional to the total area of its incident triangles. Note that the gauge is an approximation of the geodesic distance on the original surface. We run Algorithm 1 for various n < 6890 with unit step size t = 1. To initialize the method, we construct a subset X0 ⊂ X by greedily adding farthest points from the current X0 . More precisely, starting with X0 = {x1 } for some initial choice x1 , we add any element from argmax min g(x, x′ ), ′ x∈X

x ∈X0

to X0 . This step is repeated until X0 has n distinct elements, that is, X0 = {x1 , . . . , xn }. Then Algorithm 1 is initialized with πi = ξ|Vi , i ∈ [n], where (Vi )ni=1 is the Voronoi partition of X with respect to (g(·, xi ))ni=1 . We also apply the Wasserstein quantization algorithm (Euclidean k-means) for the same values of n and the same initialization π. Both algorithms converge after fewer than 100 iterations, and return a partitioning (Vi )ni=1 of the vertex set X for each n. From the partitions a 3D graph structure is reconstructed as follows. We associate each cluster Vi with its mean, that is yi :=

1 X ξ({x}) x, ξ(Vi )

i ∈ [n],

x∈Vi

and we place an edge between yi and yj , i ̸= j, whenever there exist xi ∈ Vi and xj ∈ Vj that are connected by an edge in the original graph. The resulting 3D graphs for both methods and various n are shown in Figure 1. The subject exhibits multiple (almost) self-contacts between the hands and head, as well as between both legs. We notice that the Wasserstein approach merges the regions of self-contact across all quantizations, drastically altering the topology of the induced graph. This is due to the fact that vertices of these near self-contact regions possess small Euclidean distances to vertices of the opposing contact region, making them likely to be clustered together. In contrast, the GW approach allows us to set the focus to cluster according to the geodesic (rather than Euclidean) distance, so that these self-contact regions are not merged. We deliberately chose a pose with near self-contacts for this experiment in order to illustrate the distinct behavior of the two methods. The figure also shows that the GW-based approach seems to quantize the limbs of the subject to a one-dimensional geometry faster than the Wasserstein-based approach with decreasing n. This effect can be seen well when considering

18

F. BEIER, S. ECKSTEIN

Figure 1. Progressive quantization of a 3D surface for n = 1000, 500, 300, 100, 50 together with the original. Top: GW quantizers with geodesic gauge. Bottom: Wasserstein quantizers with respect to the Euclidean metric. Insets (left): a magnified view of the subject’s right forearm and hand at n = 300, with each mesh edge colored by the cluster shared by its two endpoints – GW (green frame) versus Wasserstein (orange frame).

the results for n = 300. Here, the forearms and hands, as well as the lower legs and feet are one-dimensional lines in the GW case, but are largely still 3D in the Wasserstein case. The insets of Figure 1 make this explicit for the forearm: colouring each mesh edge by the cluster of its endpoints, the GW clusters appear as bands wrapping successive cross-sections, whereas the Wasserstein clusters are compact patches that ignore the circumference. Since a limb is thin compared to its length, the GW quantization can collapse each cross-section at almost no cost. More precisely, two points x, x′ ∈ X on opposite sides of a thin cross-section have almost identical gauge profiles, g(x, ·) ≈ g(x′ , ·). There is thus an incentive to cluster cross-sections, which roughly preserves the distances along the limb. The Wasserstein quantization objective, in contrast, sees only ambient position and forms compact three-dimensional clusters. Thereby the full thickness of a limb is retained until n is small enough that the diameter of these clusters exceeds the limb thickness. 6.2. Accelerated Pairwise Distance Computation. In this example, we show how Algorithm 1 can be used to accelerate the computation of all pairwise GW distances within a dataset of 3D shapes. We turn to the deformation-transfer dataset of Sumner and Popović [62] and extract exactly one triangulated mesh of each of the 8 classes camel, cat, elephant, horse, lion, face, head, flamingo. We pre-process each mesh and extract a gm-space similarly to Subsection 6.1. More precisely, each mesh is decimated to approximately 1000 vertices, see Figure 2. Then, for each of the 8 meshes, we extract a gm-space Xk = (Xk , gk , ξk ), where

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

19

Figure 2. A set of 8 decimated 3D surfaces from Sumner and Popović [62].

Xk ⊂ R3 are the vertex positions, gk are the normalized pairwise Dijkstra distances of vertices within the graph induced by the mesh and ξk is the normalized area-weighted vertex measure on Xk , k = 1, . . . , 8. We consider the task of computing all pairwise distances between
X1 , . . . , X8 which amounts to 82 = 28 GW computations. We use our quantization machinery to reduce the computational cost. More precisely, fixing some n < 1000, we propose to quantize each object Xk using Algorithm 1, yielding some Yk , k = 1, . . . , 8. The initialization of the algorithm is exactly done as in Subsection 6.1. Then we approximate (6.1)

GW2 (Xk , Xl ) ≈ GW2 (Yk , Yl ),

k, l = 1, . . . , 8.

Since the quantizations are much smaller than the original spaces, we expect a significant speed-up at the cost of some approximation error. Before proceeding, we remark that the latter is approximately bounded. Indeed, for arbitrary quantizers Yk , the triangle and reverse triangle inequalities yield |GW2 (Xk , Xl ) − GW2 (Yk , Yl )| ≤ |GW2 (Xk , Xl ) − GW2 (Xk , Yl )|+|GW2 (Xk , Yl ) − GW2 (Yk , Yl )| ≤ GW2 (Xk , Yk ) + GW2 (Xl , Yl ). The right-hand side is the total quantization cost of the two shapes, which is returned as a GW byproduct of Algorithm 1. It upper-bounds qGW n (Xk ) + qn (Xl ) and thus provides an a-priori 2 error bound for (6.1). As a baseline, we compute the distances (GW2 (Xk , Xl ))1≤k<l≤8 directly, which takes 178 seconds. We apply the described methodology for n = 25, 50, 75, 100, 150, 200, 300, 400. In Figure 3, we plot the total time of computing all 8 quantizations Y1 , . . . , Y8 and all their pairwise GW distances (GW2 (Yk , Yl ))1≤k<l≤8 against the mean relative error 1 28 2

X 1≤k<l≤8

|GW2 (Yk , Yl ) − GW2 (Xk , Xl )| . GW2 (Xk , Xl )

In our experiments the numerical analogue to this bound turns out to be conservative: the observed errors are usually about an order of magnitude below it.

20

F. BEIER, S. ECKSTEIN

Figure 3. Total runtime (x-axis) versus mean relative error (y-axis) for approximating all 28 pairwise GW2 distances by quantization to n points. Each marker corresponds to one value of n; the dotted red line is the direct full-resolution baseline. Table 1 GW W Constants in the two-sided bound c qW n (ξ) ≤ qn (X) ≤ C qn (ξ) for the three Euclidean gauges. Here λmin (·) and λmax (·) denote the smallest and largest eigenvalue of a matrix, Σξ is the covariance matrix of ξ, and Mξ is its second-moment matrix.

g

c

C

∥· − ·∥ ∥· − ·∥2 ⟨·, ·⟩

p −1 diam(X) p · λmin (Σξ ) 2p λmin (Σξ ) λmin (Mξ )

2 4 diam(X) p 2 λmax (Mξ )

As expected, the relative error is generally decreasing for increasing n, where the slight increase at n = 200 is likely due to the non-convexity of the GW problem, whose solver only attains local optima. In summary, already for moderate n the quantization approach reduces the runtime significantly while incurring a relative error of only around two to ten percent. 6.3. Verification of Bounds and Rates. In Section 4, under suitable conditions we have shown two-sided bounds of the form GW W c · qW n (ξ) ≤ qn (X) ≤ C · qn (ξ),

when X = (X, g, ξ) ∈ GM for X ⊂ Rd , and g ∈ {∥ · − · ∥, ∥ · − · ∥2 , ⟨·, ·⟩}, for all n ∈ N. The Wasserstein quantization value qW n (ξ) is always taken with respect to the Euclidean distance. The exact constants c, C > 0 depend on the choice of gauge g and are summarized in Table 1. To showcase the behavior of the bounds in practice, we sample N = 5 000 points x1 , . . . , xN from the uniform distribution on [0, 1]d for d = 1, 2, 3. We set X = (X, g, ξ) ∈ GM, where X ⊂ 1 PN d 2 [0, 1] is the convex hull of {x1 , . . . , xN }, g ∈ {∥ · − · ∥, ∥ · − · ∥ , ⟨·, ·⟩} and ξ = N i=1 δxi . Note

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

21

Figure 4. Verification of the two-sided bounds and the quantization rate. For each dimension d ∈ {1, 2, 3} (columns) we sample N = 5 000 points uniformly from [0, 1]d and form X = (X, g, ξ) with sampled points X, uniform measure ξ, and one of the three Euclidean gauges g (rows). Each panel shows qGW n (X) (approximated W via Algorithm 1) in red against n on log–log axes. It lies in the grey band [ c qW n (ξ), C qn (ξ) ] (constants from Table 1, qW n (ξ) approximated via Lloyd’s algorithm), and its fitted slope approximately matches the predicted rate −1/d. Here qW n (ξ) is the Wasserstein quantization value with respect to the Euclidean metric.

that X is covered by the theory of Section 4 that concerns the two-sided bounds. We proceed to run Algorithm 1 to approximately compute qGW n (X) numerically, for n = 8, 16, 32, 64, 128, 256, initialized by taking the Voronoi cells of a greedy n-furthest-point-sampling procedure on X according to the Euclidean distance, cf. Subsection 6.1. We also use the classic Lloyd algorithm, initialized with the same subsampled points, to approximate qW n (ξ). In every case, Algorithm 1 and Lloyd’s algorithm converge in under 200 iterations. Based on the W W approximations, we plot the graph of n 7→ qGW n (X), as well as the band n 7→ [cqn (ξ), Cqn (ξ)] 2 for all three gauges g ∈ {∥ · − · ∥, ∥ · − · ∥ , ⟨·, ·⟩} and all choices of dimensions d = 1, 2, 3, as log-log plots in Figure 4. In every case, the numerical approximation of the quantization value qGW n (X) lies within W (ξ)]. Note that for the scalar the numerical approximations of the interval [c qW (ξ), C q n n product with d = 1, the values lie essentially at the upper bound. This reflects Theorem 4.8, p W (ξ) 2M − qW (ξ)2 . As d = 1, M which provides the closed-form expression qGW (X) = q ξ ξ n n n W (ξ) → 0 as n → ∞, the ratio reduces to a scalar and moreover M = λ (M ). Since q max ξ ξ n p W qGW 2λmax (Mξ ) = C, so the upper bound is asymptotically attained. n (X)/qn (ξ) converges to

22

F. BEIER, S. ECKSTEIN

We use a large sample size N = 5 000 ≫ n so that qGW n (X) serves as a proxy for the quantization value of the continuous uniform measure on [0, 1]d . Accordingly, the log–log plots display an approximately linear decay, and the fitted slopes match the predicted rate −1/d to within ≈ 0.02 across all gauges and dimensions. This reflects the asymptotic rate −1/d established in Corollaries 4.4 and 4.6. qGW n (X) ≍ n 6.4. Pruning a Neural Network. Pruning neural networks is a way of reducing the size (and thus the computational burden) of a large neural network while trying to retain a high performance. In this example, we show how a feedforward neural network can be encoded as a gm-space and how the respective quantization leads to a smaller network which performs better than baseline pruning procedures. To this end, we work with a simple three-hiddenlayer neural network of the form f (x) = W4 ◦ φ ◦ W3 ◦ φ ◦ W2 ◦ φ ◦ W1 (x), where x ∈ Rd is the input to the neural network, W1 , W2 , W3 , W4 are linear maps (i.e., matrices), and φ is an activation function applied element-wise. Most pruning methods work layer by layer. A simple pruning baseline is to delete neurons based on certain characteristics, usually related to the magnitudes of the attached weights (see, e.g. [25, 32, 36]) or activation vectors (see, e.g., [63]). Beyond just deletion of neurons, pruning can be based on clustering each layer’s neurons and suitably combining the weights of neurons which share a cluster (see, e.g., [1, 44, 58]). While pruning methods are sometimes combined with retraining, in this example we focus purely on data-free methods, and hence neither retraining nor activation vectors are used. In the following, we show how to encode and prune a neural network as a whole using the presented GW machinery. To this end, we define the gm-space X = ([N ], g, ξ), where N is the total number of neurons (including input and output neurons) in the network and ξ is the uniform distribution. In the simple experiment (where the MNIST dataset is used which consists of images of handwritten digits), we have N = 784 + 3 · h + 10, where 784 is the input dimension (28 × 28 pixel images), h ∈ {50, 100} is the hidden dimension for our experiments, and 10 is the output dimension, corresponding to 10 digits. The gauge g is built via the weight matrices of the network as follows: For l, ˜l ∈ {0, 1, 2, 3, 4} with l < ˜l, set W (l,l̃) := Wl̃ Wl̃−1 . . . Wl+1 and W (l̃,l) := (W (l,l̃) )⊤ . Then, if i is the a-th neuron in layer l and j (l,l̃) is the b-th neuron in layer ˜l, where l < ˜l, we set g(i, j) := λ(|l−l̃|−1)/2 W , where λ ∈ (0, 1) is b,a

a factor which we introduce to increase the relative importance of short-range interactions in the network. In the following, we use λ = 0.3, but we observed similar results for λ ∈ [0.1, 0.5]. If l = ˜l and i ̸= j, we set g(i, j) = 0. Regarding the diagonal g(i, i), we make the following observations: First, in a clustering sense, each input and each output neuron is unique, and hence we set their gauge value g(i, i) each to a unique value separated sufficiently from all other diagonal values. This ensures that no input or output neuron should ever be clustered together with any other neuron. Similarly, so that no neurons are matched across different hidden layers, each hidden layer has a unique value for their corresponding neurons on the diagonal, each sufficiently separated from all other values. Effectively, what these diagonal values as a whole are doing is to restrict the forms of the support which the measures πi in Algorithm 1 can take. Finally, we symmetrize the gauge by setting g(j, i) := g(i, j).

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

23

Figure 5. Comparison of different ways to reduce the size of the neural network (three-hidden-layer feedforward neural network trained on MNIST) by reducing each hidden dimension to a width multiplier times the initial width. The baselines are pruning (gray lines, arising from deleting those neurons per layer whose incoming weights have lowest L1 norm) and layer-wise clustering of neurons (blue lines, cf. [44]). The method introduced in this paper is shown as the green lines, where the network as a whole is modeled as a gm-space and quantized with Algorithm 1. The confidence bands are over 10 independent random seeds.

To be precise, for neurons i ∈ [N ], we set l(i) ∈ {0, 1, 2, 3, 4} as its layer, and we set p(i) ∈ {1, 2, 3, . . .} as the neuron’s position within the layer. For i, j ∈ [N ], the gauge is given by   vi i = j and l(i) ∈ {0, 4}    y i = j and l(i) ∈ {1, 2, 3} l(i) g(i, j) :=  0 l(i) = l(j) and i ̸= j    (l(i),l(j)) λ(|l(i)−l(j)|−1)/2 W l(i) ̸= l(j) p(j),p(i)

where vi , yl(i) ∈ R are just sufficiently different values on the diagonal. We run Algorithm 1 to quantize X. To this end, we note that we merely reduce the size of the hidden layers, so we prune the model size to n = 784 + 3 · ⌊α · h⌋ + 10, where α ∈ [0, 1] is a multiplier for the width of hidden layers. To initialize Algorithm 1, we use the layer-by-layer

24

F. BEIER, S. ECKSTEIN

hierarchical clustering of neurons as given by [44]. We then run Algorithm 1 with step size 1, but note that the results are practically the same when using a line search instead. To ensure a fair comparison with the baseline methods, we want to avoid empty clusters in Algorithm 1, for which we use standard cluster splitting methods, cf. [64, Section 8.2.2]. We denote the output of Algorithm 1 by ([n], G, υ). Due to the aforementioned modeling choice of the gauge (effectively restricting πi ), we can still assign each neuron i ∈ [n] uniquely to a layer. To extract the pruned neural network from this output, we first fit a gauge of neural-networkform (as g, but for a smaller network) to G in L2 (υ ⊗ υ), which we denote by G̃.3 Finally, for two neurons i, j ∈ [n] in layers l and l + 1, the corresponding entry of the weight matrix W̃l+1 υi of the pruned network is given by ξ({1}) G̃(i, j).4 The results for networks with varying activation functions and sizes are shown in Figure 5, which shows results for networks trained on MNIST. We see that the introduced method using GW quantization performs best throughout all four regimes (different sizes and activation functions), particularly leading to better performance at more aggressive compression. While this example is certainly only a first step in this direction (e.g., no retraining is used and the modeling is agnostic to the activation function), it showcases the potential of the modeling freedom provided by GW to encode neural networks and the relevance of quantization in such settings. Appendix A. Auxiliary Lemmas. We provide two auxiliary lemmas which are known, but stated here for clear reference in our setting. Lemma A.1. Let ξ ∈ P(X) and (ci )ni=1 with non-negative c1 , . . . , cn ∈ L1 (ξ). Then n Z X

min

π1 ,...,π Pnn ∈M+ (X) i=1 i=1 πi =ξ

Z ci dπi (x) =

X

min cj dξ.

X j=1,...,n

Pn

A collection of measures (πi )ni=1 ⊂ M+ (X) with if πi is concentrated on

i=1 πi = ξ attains the minimum if and only

Ai := {x ∈ X : ci (x) = min cj (x)}, j=1,...,n

for all i = 1, . . . , n.

In particular, if (Vi )ni=1 is a Voronoi partition of X according to (ci )ni=1 , (ξ|Vi )ni=1 is a solution and it is unique if and only if ξ(Ai ∩ Aj ) = 0, i ̸= j. Proof. We abbreviate c := minj=1,...,n cjP . Since ci ∈ L1 (ξ) and 0 ≤ c ≤ ci , we also have 1 n c ∈ L (ξ). Now let (πi )i=1 ⊂ M+ (X) with ni=1 πi = ξ be arbitrary. Then (A.1)

n Z X i=1

3

X

ci dπi ≥

n Z X i=1

X

Z c dπi =

c dξ. X

We perform this L2 -minimization via the Adam optimizer, which appears to converge quickly and robustly. υi The inclusion of the term ξ({1}) is important to retain the scale of the weights in each layer. It encodes how many initial neurons are summarized by the new node i, which is not included in the gauge value (which is simply an average across initial gauge values). Of course, we could replace ξ({1}) by 1/N . 4

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

25

It holds ci − c ≥ 0, and X \ Ai = {x : ci (x) − c(x) > 0}. Thus, equality in (A.1) holds if and only if πi (X \Ai ) = 0, i.e. πi is concentrated on Ai , for every i = 1, . . . , n. By construction, the Voronoi partition (Vi )ni=1 is a disjoint partition of X with Vi ⊂ Ai , so in particular (ξ|Vi )ni=1 is feasible and attains the lower bound in (A.1). We turn to the uniqueness claim. First, assume ξ(Ai ∩ Aj ) = 0 for all i ̸= j and let (πi )ni=1 be any solution. It holds (A.2)

ξ(B ∩ Vi ) =

n X

πj (B ∩ Vi ) = πi (B ∩ Vi ) +

j=1

X

πj (B ∩ Vi ),

B ⊂ X measurable.

j̸=i

For j ̸= i, we have Vi ∩ Aj ⊂ Ai ∩ Aj , so that πj (Vi ) = πj (Vi ∩ Aj ) ≤ πj (Ai ∩ Aj ) ≤ ξ(Ai ∩ Aj ) = 0. Hence, together with (A.2), we obtain πi |Vi = ξ|Vi . In the following, we establish that πi is concentrated on Vi . If x ∈ Ai \ Vi , then x ∈ Vj ⊂ Aj for some j ̸= i and thus x ∈ Ai ∩ Aj . Hence, [ Ai \ V i ⊂ (Ai ∩ Aj ), j̸=i

which is a ξ-null set. As πi ≤ ξ and πi is concentrated on Ai = Vi ∪ (Ai \ Vi ), it follows that πi is concentrated on Vi and πi = ξ|Vi , for all i = 1, . . . , n. Conversely, assume ξ(Ak ∩ Al ) > 0 for some fixed k ̸= l. Since (Vi )ni=1 is a partition of X, there exists i ∈ {1, . . . , n} such that M := Ak ∩ Al ∩ Vi satisfies ξ(M ) > 0. Choose j ∈ {k, l} \ {i} and set ′ := ξ|Vm , m ∈ {1, . . . , n} \ {i, j}. πm Pn ′ ′ ′ )n Since M ⊂ Vi , we obtain (πm m=1 ⊂ M+ (X) with m=1 πm = ξ. Each πm is concentrated ′ )n on Am , since Vm ⊂ Am and M ⊂ Aj as well as M ⊂ Vi ⊂ Ai . By the first part, (πm m=1 is a n solution, and it differs from (ξ|Vm )m=1 as ξ(M ) > 0, as desired.

πi′ := ξ|Vi − 21 ξ|M ,

πj′ := ξ|Vj + 12 ξ|M ,

The following is a well-known statement regarding the arithmetic mean as the L2 barycenter. Lemma A.2. Let A be a Polish space, α R∈ M+ (A) with α(A) > 0 and f ∈ L2 (A, α; Rd ). 1 d Then, for its associated L2 -mean m := α(A) A f dα ∈ R , we have ∥f − m∥L2 (α) = inf ∥f − c∥L2 (α) ,

(A.3)

c∈Rd

and m is the unique minimizer. Moreover, it holds ∥f − a0 ∥2L2 (α) = ∥f − m∥2L2 (α) + ∥m − a0 ∥2 α(A), Proof. For any a0 ∈ Rd , expanding the square gives Z 2 ∥f − a0 ∥L2 (A,α) = ∥f (a) − a0 ∥2 dα(a) A

for any a0 ∈ Rd .

26

F. BEIER, S. ECKSTEIN

Z 

 ∥f (a) − m∥2 + 2⟨f (a) − m, m − a0 ⟩ + ∥m − a0 ∥2 dα(a) A E DZ (f (a) − m) dα(a), m − a0 + ∥m − a0 ∥2 α(A). = ∥f − m∥2L2 (A,α) + 2

=

A

R

R By definition of m, it holds A (f (a) − m) dα(a) = A f dα − m α(A) = 0, so the middle term vanishes and the asserted identity follows. It holds ∥m − a0 ∥2 α(A) ≥ 0 and equality if and only if a0 = m, this also yields (A.3) and the uniqueness of m. Appendix B. Proofs for Section 3. B.1. Fully Supported Representatives in GMn . Proof of Proposition 3.1. Let Y = ([n], G, υ) ∈ GMn and set Y1 = {i ∈ [n] : υi > 0}. We assume m := |Y1 | < n, otherwise we are done. Firstly, Y1 := (Y1 , G|Y1 ×Y1 , υ|Y1 ) is homomorphic to Y by construction. Moreover, by renaming Y1 according to any bijection T1 : Y1 → [m], Y1 is further homomorphic to some Y2 = ([m], G(2) , υ (2) ) with υ (2) ∈ P([m]) fully supported. It suffices to find a to-Y2 homomorphic gm-space of the form Ỹ = ([m + 1], G̃, υ̃), with fully supported υ̃ ∈ P([m + 1]). From there the result follows by induction. For this let b = (2) (2) (G1,m , . . . , Gm,m )⊤ ∈ Rm and set G̃ =

G(2) b (2) ⊤ b Gm,m

! (m+1)×(m+1)

∈R

 and υ̃ =

(2) (2)  υm υm ⊤ (2) (2) υ1 , . . . , υm−1 , , .

2

2

By construction, the map T : [m + 1] → [m] given by T |[m] = Id[m] , and T (m + 1) = m, defines a homomorphism from Ỹ to Y2 . Hence, GW2 (Y2 , Ỹ) = 0, so that Y is homomorphic to Ỹ. Repeating the construction inductively for m + 2, . . . , n, we see that Y is homomorphic to a space defined on [n] with fully supported measure. B.2. Optimality of the Block-Mean Gauge. Proof of Theorem 3.2. Consider a fixed π ∈ P(X × [n]) with (PX )# π = ξ. Clearly, Ĝ is symmetric. Furthermore, it holds inf QGW (G, (πi )ni=1 ) = inf G∈Rn×n G∈Rn×n =

(B.1)

=

n X i,i′ =1 n X i,i′ =1

n ZZ X i,i′ =1

|g(x, x′ ) − Gi,i′ |2 dπi (x) dπi′ (x′ )

X2

inf ∥g − Gi,i′ ∥2L2 (πi ⊗π ′ )

Gi,i′ ∈R

i

∥g − Ĝi,i′ ∥2L2 (πi ⊗π ′ ) , i

where the last equality follows by Lemma A.2, since Ĝi,i′ is exactly the L2 mean of g with respect to πi ⊗ πi′ , i, i′ ∈ [n]. Since the unconstrained minimizer is symmetric, it also solves the problem when constrained to symmetric matrices. Consider the case π ∈ Πo (X, Y) for the

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

27

given Y = ([n], G, υ). Then GW22 (X, Y) = QGW (G, (πi )ni=1 ) =

n X

∥g − Gi,i′ ∥2L2 (πi ⊗π ′ ) i

i,i′ =1

 n  X 2 2 = ∥g − Ĝi,i′ ∥L2 (πi ⊗π ′ ) + πi (X)πi′ (X)(Gi,i′ − Ĝi,i′ ) i {z } | i,i′ =1 =υi υi′ ZZ |g(x, x′ ) − Ĝi,i′ |2 dπ(x, i) dπ(x′ , i′ ) + ∥G − Ĝ∥2L2 (υ⊗υ) =

(X×[n])2 ≥ GW22 (X, Ŷ) + ∥G − Ĝ∥2L2 (υ⊗υ)

where the equality in the second line follows again from Lemma A.2 and the final estimate is due to π ∈ Π(ξ, υ). B.3. Optimal Couplings Concentrate on Voronoi Cells. First, we show that QGW (G, ·) attains its minimum over Π(ξ, ∗n ), for which we require an intermediary result which we present in a slightly generalized form that is required later. Lemma B.1. Let X be a Polish space, ξ ∈ P(X) and consider two sequences (π (k) )k∈N , (σ (k) )k∈N ⊂ P(X × [n]) with (PX )# π (k) = (PX )# σ (k) = ξ, k ∈ N, such that π (k) ⇀ π and σ (k) ⇀ σ. Then, for every h ∈ L1 (ξ ⊗ ξ) and all i, i′ ∈ [n], ZZ ZZ (k) (k) ′ k→∞ ′ h(x, x′ ) dπi (x) dσi′ (x′ ). h(x, x ) dπi (x) dσi′ (x ) −−−→ X2

X2

(k)

(k)

Proof. Since πi (A) = π (k) (A × {i}) ≤ ξ(A) and, likewise, σi′ (A) ≤ ξ(A) for all measur(k) (k) able A ⊂ X, it holds πi ⊗ σi′ ≤ ξ ⊗ ξ. Thus ZZ

(k)

(B.2) X2

f dπi

(k)

dσi′

≤ ∥f ∥L1 (ξ⊗ξ)

for all f ∈ L1 (ξ ⊗ ξ), k ∈ N.

As the marginal projection is weakly continuous, the limits satisfy (PX )# π = (PX )# σ = ξ, so (k) (k) πi ⊗ σi′ ≤ ξ ⊗ ξ and (B.2) also holds with πi , σi′ in place of πi , σi′ . First, let h̃ ∈ Cb (X × X), the function ((x, j), (x′ , j ′ )) 7→ h̃(x, x′ )1{i} (j)1{i′ } (j ′ ) is bounded and continuous on (X × [n])2 . Hence, due to π (k) ⊗ σ (k) ⇀ π ⊗ σ see [12, Prop. 2.7.8], we obtain ZZ ZZ (k) (k) k→∞ (B.3) h̃ dπi dσi′ −−−→ h̃ dπi dσi′ . X2

X2

Now, let h ∈ L1 (ξ ⊗ ξ) and ε > 0. As Cb (X × X) is dense in L1 (ξ ⊗ ξ), fix h̃ ∈ Cb (X × X) with ∥h − h̃∥L1 (ξ⊗ξ) < ε. Applying (B.2) to f = h − h̃, for both π (k) and π, ZZ X2

(k) (k) h dπi dσi′ −

ZZ

ZZ h dπi dσ

X2

i′

≤ X2

(k) (k) (h − h̃) dπi dσi′

ZZ (h̃ − h) dπi dσi′

+ X2

28

F. BEIER, S. ECKSTEIN

ZZ

ZZ (k) (k) h̃ dπi dσi′ − h̃ dπi dσi′ 2 2 X X ZZ ZZ (k) (k) ≤ 2ε + h̃ dπi dσi′ − h̃ dπi dσi′ . +

X2

Taking lim supk→∞ and using (B.3) yields ZZ ZZ (k) (k) lim sup h dπi dσi′ − k→∞

X2

X2

h dπi dσi′ ≤ 2ε. X2

Since ε > 0 was arbitrary, the left-hand side is zero, as desired. (n)

Proposition B.2. Let X = (X, g, ξ) ∈ GM and G ∈ Gsym . Then QGW (G, ·) attains its minimum over Π(ξ, ∗n ). Proof. The set Π(ξ, ∗n ) is weakly compact [70, Cor. 5.21]. Writing n ZZ X

QGW (G, (πi )ni=1 ) =

(g(x, x′ ) − Gi,i′ )2 dπi (x) dπi′ (x′ )

X2

i,i′ =1

with (g − Gi,i′ )2 ∈ L1 (ξ ⊗ ξ), Lemma B.1 shows that QGW (G, ·) is weakly continuous on Π(ξ, ∗n ). The claim follows by the Weierstrass theorem. Now that we established existence of minimizers, we want to characterize them. Unfortunately, the objective (πi )ni=1 7→ QGW (G, (πi )ni=1 ) is quadratic and its minimization does not admit closed form as in the linear Wasserstein setting. Following the GW literature [48, 56, 68], a common approach to handle such problems, is via their bilinear relaxation. More precisely, we define Q̃GW (G, (πi )ni=1 , (π̃i )ni=1 ) :=

n ZZ X i,i′ =1

(g(x, x′ ) − Gi,i′ )2 dπi (x) dπ̃i′ (x′ ),

X2

which is linear in (πi )ni=1 and (π̃i )ni=1 , respectively. For some fixed π̃ ∈ P(X × [n]) with (PX )# π̃ = ξ, we define (π̃) ci (x) :=

n Z X i′ =1

(g(x, x′ ) − Gi,i′ )2 dπ̃i′ (x′ ), (π̃)

Since π̃i ≤ ξ and g ∈ L2 (ξ ⊗ ξ), it holds ci (B.4)

i ∈ [n].

X

∈ L1 (ξ). Using this new notation, we have

Q̃GW (G, (πi )ni=1 , (π̃i )ni=1 ) =

n Z X i=1

X

(π̃)

ci (x) dπi (x).

In our setting, the rewriting (B.4) casts the linear minimization as a linear assignment problem, whose solutions the following proposition characterizes via Voronoi partitions.

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

29 (n)

n

a fixed gauge matrix G ∈ Gsym and (π̃i )i=1 ⊂ M+ (X) with Pn Proposition B.3. Consider n π̃ = ξ. Then, (π ) minimizes i i i=1 i=1 inf

(πi P )n i=1 ⊂M+ (X) n i=1 πi =ξ

Q̃GW (G, (πi )ni=1 , (π̃i )ni=1 )

if and only if πi is concentrated on (π̃)

(π̃)

Ai := {x ∈ X : ci (x) = min cj (x)}, j

for every i ∈ [n]. (π̃)

In particular, (ξ|Vi )ni=1 is a solution for a Voronoi partition (Vi )ni=1 of X according to (ci )ni=1 . It is the unique solution if and only if ξ(Ai ∩ Aj ) = 0 for all i ̸= j. Proof. This directly follows from (B.4) and Lemma A.1. Going via the bilinear relaxation allows us to weakly characterize the solutions of the original quadratic problem. The idea is to show that π̂, as a minimizer of the quadratic problem, also minimizes the linear function Q̃GW (G, ·, (π̂i )ni=1 ) from where the result follows by Proposition B.3. Proof of Theorem 3.3. Note that existence is covered by Proposition B.2. For brevity, we use the notation B(π, π̃) := Q̃GW (G, (πi )ni=1 , (π̃i )ni=1 ). Let (Vi )ni=1 be the partition from Proposition B.3 and let γ ∈ P(X × [n]) be defined by γ(· × {i}) = ξ|Vi . According to the proposition, γ is a solution to (B.5)

argmin

B(π, π̂).

(πi P )n i=1 ⊂M+ (X) n i=1 πi =ξ

Hence, δ := B(γ, π̂) − B(π̂, π̂) ≤ 0. We show δ = 0. For t ∈ (0, 1), we set π (t) = π̂ + t(γ − π̂). Clearly, (PX )# π (t) = ξ for all t ∈ (0, 1). Due to the bilinearity and symmetry of B, we obtain B(π (t) , π (t) ) = B(π̂, π̂) + 2t

B(π̂, γ − π̂) | {z }

+t2 B(γ − π̂, γ − π̂)

=B(π̂,γ)−B(π̂,π̂)=δ

Minimality of π̂ ensures B(π (t) , π (t) ) ≥ B(π̂, π̂), so that 0 ≤ 2δ + tB(γ − π̂, γ − π̂). As δ ≤ 0, taking the limit t → 0 yields 0 = δ = B(γ, π̂) − B(π̂, π̂). Hence, π̂ is a minimizer of (B.5), so that Proposition B.3 yields that each π̂i is concentrated on Ai . If additionally ξ(Ai ∩ Aj ) = 0 for all i ̸= j, the uniqueness part of Proposition B.3 gives π̂i = ξ|Vi , i ∈ [n], as desired. B.4. Existence and Structure of GW Quantizers. Proof of Theorem 3.4. Let X = (X, g, ξ) ∈ GM be arbitrary. We first show existence. Let (k) = ([n], G(k) , υ (k) ), k ∈ N, whose k∈N ⊂ GMn be a minimizing sequence of (3.2) with Y

(Y(k) )

30

F. BEIER, S. ECKSTEIN

measures are fully supported, see Proposition 3.1, and let π (k) ∈ Πo (X, Y(k) ), k ∈ N. Define Ŷ(k) := ([n], Ĝ(k) , υ (k) ) ∈ GMn with ZZ 1 (k) (k) (k) g(x, x′ ) dπi (x) dπi′ (x′ ), i, i′ ∈ [n]. Ĝi,i′ := (k) (k) X2 υi υi′ By Theorem 3.2, it holds (k)

(k)

2 GW22 (X, Y(k) ) = QGW (G(k) , (πi )ni=1 ) ≥ QGW (Ĝ(k) , (πi )ni=1 ) ≥ GW22 (X, Ŷ(k) ) ≥ qGW n (X) ,

so that (k)

2 (k) n qGW n (X) = lim QGW (Ĝ , (πi )i=1 ) k→∞

and (Ŷ(k) )k∈N is also a minimizing sequence of (3.2). Since (π (k) )k∈N is a sequence in the weakly compact set Π(ξ, ∗n ), we can select a converging subsequence, again denoted by (π (k) )k . Let π ∈ P(X × [n]) be its limit and set υ := (P[n] )# π which is the weak limit of the corresponding subsequence of marginals (υ (k) )k∈N . We set I := {i ∈ [n] : υi > 0} and note that πi (X) = υi = 0, i.e., πi = 0, for i ∈ / I. For i, i′ ∈ I, Lemma B.1 applied with h = g yields ZZ 1 (k) k→∞ Ĝi,i′ −−−→ Ĝi,i′ := g(x, x′ ) dπi (x) dπi′ (x′ ). υi υi′ X2 (k)

(k)

(k)

In the context of Lemma A.2 with d = 1, f = g and α = πi ⊗ πi′ we get m = Ĝi,i′ , so applying the lemma with a0 = 0 yields ZZ ZZ (k)
2 (k) (k) (k) (k) (k) (k) (k)
2 g 2 dπi dπi′ − Ĝi,i′ υi υi′ g − Ĝi,i′ dπi dπi′ = 2 X2 Z ZX
2 k→∞ −−−→ g 2 dπi dπi′ − Ĝi,i′ υi υi′ , i, i′ ∈ I, X2

where the convergence of the first term holds again by Lemma B.1. Thus, we arrive at ZZ X (k) n (k)
2 (k) (k) GW 2 (k) qn (X) = lim QGW (Ĝ , (πi )i=1 ) ≥ lim g − Ĝi,i′ dπi dπi′ k→∞

=

i,i′ ∈I

X Z Z i,i′ ∈I

X2

2

2

k→∞

g dπi dπi′ − Ĝi,i′ υi υi′

X2

 =

X ZZ i,i′ ∈I

g − Ĝi,i′

2

dπi dπi′ ,

X2

where the last equality follows as before, from Lemma A.2. Now set Ŷ := ([n], Ĝ, υ) ∈ GMn , where Ĝ is extended by zero outside of I × I. As πi = 0 for i ∈ / I, the right-hand side above equals the cost of the coupling π ∈ Π(ξ, υ) between X and Ŷ, so that ZZ
2 GW 2 2 qn (X) ≥ g(x, x′ ) − Ĝi,i′ dπ(x, i) dπ(x′ , i′ ) ≥ GW22 (X, Ŷ) ≥ qGW n (X) , (X×[n])2

where the last inequality holds since Ŷ ∈ GMn . Hence, all inequalities are equalities and Ŷ is a solution of (3.2).

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

31

We turn to the characterization. Let Ŷ = ([n], Ĝ, υ̂) be an arbitrary solution of (3.2), with fully supported υ̂ ∈ P([n]). Let X = (X, g, ξ) ∈ GM and let π̂ ∈ Πo (X, Ŷ). Note that π̂i (X) = υ̂i > 0 for all i ∈ [n]. (n) We show Item (i). Define G̃ ∈ Gsym by ZZ 1 G̃i,i′ := g(x, x′ ) dπ̂i (x) dπ̂i′ (x′ ), i, i′ ∈ [n], υ̂i υ̂i′ 2 X and set Ỹ := ([n], G̃, υ̂) ∈ GMn . Since π̂ ∈ Πo (X, Ŷ), Theorem 3.2 yields GW22 (X, Ŷ) ≥ GW22 (X, Ỹ) + ∥Ĝ − G̃∥2L2 (υ̂⊗υ̂) . 2 2 2 As Ŷ solves (3.2), it holds GW22 (X, Ỹ) ≥ qGW n (X) = GW2 (X, Ŷ), so that ∥Ĝ − G̃∥L2 (υ̂⊗υ̂) = 0.

Since υ̂ is fully supported, this gives Ĝ = G̃, as desired. We show Item (ii). With the gauge matrix Ĝ of Ŷ fixed, since π̂ ∈ Πo (X, Ŷ), it holds 2 QGW (Ĝ, (π̂i )ni=1 ) = GW22 (X, Ŷ) = qGW n (X) .

Thus, (π̂i )ni=1 is a solution to inf

(πi P )n i=1 ⊂M+ (X) n i=1 πi =ξ

QGW (Ĝ, (πi )ni=1 ).

(π̂)

Hence, Theorem 3.3 yields that π̂i is concentrated on Ai for every i ∈ [n], where ĉi = ci taken with respect to the gauge matrix Ĝ, which is Item (ii). We show Item (iii). Since π̂ ∈ Πo (X, Ŷ), it holds n ZZ X

GW22 (X, Ŷ) =

g − Ĝi,i′

2

is

dπ̂i dπ̂i′ .

X2

i,i′ =1

By Item (i), Ĝi,i′ is the mean of g with respect to the measure π̂i ⊗ π̂i′ , so that Lemma A.2 applied with a0 = 0 yields ZZ ZZ
2 g − Ĝi,i′ dπ̂i dπ̂i′ = g 2 dπ̂i dπ̂i′ − Ĝ2i,i′ υ̂i υ̂i′ , i, i′ ∈ [n]. X2

X2

Pn

Summing over i, i′ ∈ [n] and using GW22 (X, Ŷ) =

i=1 π̂i = ξ, we obtain

ZZ

2

′

′

g (x, x ) dξ(x) dξ(x ) − X2

n X

Ĝ2i,i′ υ̂i υ̂i′ ,

i,i′ =1

which is Item (iii). Finally, assume additionally that ξ(Ai ∩ Aj ) = 0 for all i ̸= j. As shown in Item P (ii), the family (π̂i )ni=1 is a solution to inf QGW (Ĝ, (πi )ni=1 ) over all (πi )ni=1 ⊂ M+ (X) with ni=1 πi = ξ, so that the second part of Theorem 3.3 yields π̂i = ξ|Vi , i ∈ [n], for the Voronoi partition (Vi )ni=1

32

F. BEIER, S. ECKSTEIN

of X with respect to (ĉi )ni=1 . Since (Vi )ni=1 is a disjoint partition of X, the map T : X → [n], T (x) = i if x ∈ Vi , is well-defined and it holds π̂ =

n X

ξ|Vi ⊗ δi = (Id, T )# ξ.

i=1

P In particular, υ̂i = π̂i (X) = ξ(Vi ) for all i ∈ [n], i.e., υ̂ = ni=1 ξ(Vi ) · δi . Inserting π̂i = ξ|Vi into Item (i) gives ZZ 1 Ĝi,i′ = g(x, x′ ) dξ(x) dξ(x′ ), i, i′ ∈ [n], ξ(Vi )ξ(Vi′ ) Vi ×Vi′ as desired. Appendix C. Proofs for Section 4. C.1. An Upper Bound for Lipschitz Functions of Metrics. Proof of Proposition 4.1. Let ε > 0 be arbitrary. By (2.3), choose Y = {y1 , . . . , yn } ⊂ X with Z 2 qW (ξ) ≥ d2X (x, Y ) dξ(x) − ε, n X

and let h := φ ◦ dX |Y ×Y . In addition, let υ = T# ξ ∈ P(Y ) be defined by T : X → Y, x 7→ yi ,

where i = min{j : dX (x, yj ) ≤ dX (x, yk ) for all k ∈ [n]}.

Let Y := (Y, h, υ) and note that qGW n (X) ≤ GW2 (X, Y) since Y ∈ GMn . We apply the triangle inequality to obtain dX (x, x′ ) − dX (T (x), T (x′ ))
≤ dX (x, T (x)) + dX (T (x), T (x′ )) + dX (T (x′ ), x′ ) − dX (T (x), T (x′ )) =dX (x, T (x)) + dX (x′ , T (x′ )), for any x, x′ ∈ X. An analogous application also yields the same estimate for dX (T (x), T (x′ ))− dX (x, x′ ), so that in total (C.1)

|dX (x, x′ ) − dX (T (x), T (x′ ))| ≤ dX (x, T (x)) + dX (x′ , T (x′ )).

For the coupling π := (Id, T )# ξ ∈ Π(ξ, υ), we obtain ZZ 2 GW2 (X, Y) ≤ (φ(dX (x, x′ )) − φ(dX (y, y ′ )))2 dπ(x, y) dπ(x′ , y ′ ) (X×Y )2 ZZ 2 ≤L (dX (x, x′ ) − dX (y, y ′ ))2 dπ(x, y) dπ(x′ , y ′ ) 2 (X×Y ) ZZ = L2 (dX (x, x′ ) − dX (T (x), T (x′ )))2 dξ(x) dξ(x′ ). X2

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

33

Now applying (C.1) together with the elementary inequality (a + b)2 ≤ 2(a2 + b2 ) gives ZZ GW 2 2 2 (dX (x, T (x)) + dX (x′ , T (x′ )))2 dξ(x) dξ(x′ ) qn (X) ≤ GW2 (X, Y) ≤ L 2 X ZZ 2 d2X (x, T (x)) + d2X (x′ , T (x′ )) dξ(x) dξ(x′ ) ≤ 2L 2 X Z 2 d2X (x, T (x)) dξ(x) ≤ 4L X 2 2 ≤ 4L2 qW n (ξ) + 4L ε. W Letting ε → 0 and taking square roots gives qGW n (X) ≤ 2L qn (ξ), as desired.

C.2. Lower Bounds for the Euclidean and Squared-Euclidean Gauge. For the proof of Theorem 4.3, we require two auxiliary lemmas, and the following additional notation. For a measure µ ∈ P(Rd ) and f : Rd → R, we write Z ZZ 2
2 1 f dµ 2 = (C.2) Varµ (f ) := f − f (x) − f (x′ ) dµ(x) dµ(x′ ), 2 Rd ×Rd L (µ) Rd R 2 where the second equality follows by expanding the squares (both sides equal f dµ − R 2 ( f dµ) ). Lemma C.1. Let µ, ν ∈ P(Rd ) and f ∈ L2 (µ ⊗ ν). Then for every c ∈ R, Z
2 ∥f − c∥L2 (µ⊗ν) ≥ Varµ f (·, y) dν(y). Rd

Proof. Let c ∈ R be arbitrary. By Lemma A.2, ∥f (·, y) − c∥2L2 (µ) ≥ Varµ f (·, y) for ν-a.e. fixed y. Integrating over y with respect to ν yields the claim. Lemma C.2. Let µ, ν ∈ P(Rd ), where µ has finite fourth and ν finite second moment, with covariance matrices Σµ , Σν . Then Z
Varµ ∥· − y∥2 dν(y) ≥ 4 tr(Σµ Σν ). Rd

R R Proof. Set m2 := Rd ∥x∥2 dµ(x), m(µ) := Rd x dµ(x) and cµ := Rd ∥x∥2 − m2 x dµ(x). Expanding ∥x − y∥2 = ∥x∥2 − 2⟨x, y⟩ + ∥y∥2 and using that Varµ is invariant under adding the (in x) constant ∥y∥2 ,

q(y) := Varµ ∥· − y∥2 = Varµ ∥·∥2 − 2⟨·, y⟩ . R

The function x 7→ ∥x∥2 − 2⟨x, y⟩ has µ-mean m2 − 2⟨m(µ), y⟩, so by the definition of the variance, Z  2
q(y) = ∥x∥2 − m2 − 2⟨x − m(µ), y⟩ dµ(x). Rd

R We expand the square. Firstly, due to Rd ∥x∥2 − m2 dµ(x) = 0, it holds Z
∥x∥2 − m2 ⟨x − m(µ), y⟩ dµ(x) = ⟨cµ , y⟩. Rd

34

F. BEIER, S. ECKSTEIN

Secondly, we get

R

2 ⊤ Rd ⟨x − m(µ), y⟩ dµ(x) = y Σµ y. Thus, we obtain

q(y) = Varµ ∥·∥2 − 4⟨cµ , y⟩ + 4 y ⊤ Σµ y. Thus, q is a polynomial of degree (at most) two in y with Hessian ∇2 q = 8Σµ . For polynomials of degree two, R the second-order Taylor expansion is exact, so that expanding q around m(ν) and using Rd (y − m(ν)) dν(y) = 0 yields (cf. [55, Theorem 1.5]), Z Z
1 (y − m(ν))⊤ ∇2 q (y − m(ν)) dν(y) q dν = q(m(ν)) + 2 Rd Rd Z

1 = q(m(ν)) + tr ∇2 q (y − m(ν))(y − m(ν))⊤ dν(y) 2 Rd = q(m(ν)) + 4 tr(Σµ Σν ).
Since q(m(ν)) = Varµ ∥· − m(ν)∥2 ≥ 0 as a pointwise variance, the claim follows. We finally turn to the proof of Theorem 4.3. Proof of Theorem 4.3. Let Y = ([n], G, υ) ∈ GMn be arbitrary with fully supported υ ∈ P([n]),R see Proposition 3.1, and let π ∈ Πo (X,R Y). We set π̄i := πi /υi ∈ P(X), with mean mi := X x dπ̄i (x) and covariance matrix Σi := X (x − mi )(x − mi )⊤ dπ̄i (x), i ∈ [n]. Note that by closed convexity of X, mi ∈ X. Since πi ≤ ξ, all appearing moments are finite. Due to π ∈ Πo (X, Y), it holds GW22 (X, Y) =

n X

2

υi υi′ g − Gi,i′ L2 (π̄ ⊗π̄ ′ ) . i

i,i′ =1

i

Applying Lemma C.1 to each summand with f = g, µ = π̄i , ν = π̄i′ and c = Gi,i′ , we obtain Z Z n n X X

2 (C.3) GW2 (X, Y) ≥ υi υi′ Varπ̄i g(·, y) dπ̄i′ (y) = υi Varπ̄i g(·, y) dξ(y), i,i′ =1

X

i=1

X

P where, in the last equality, we used ni′ =1 υi′ π̄i′ = ξ. We show Item (i), i.e. let g = ∥· − ·∥2 . By Lemma C.2 applied with µ = π̄i and ν = ξ, and since tr(AB) ≥ λmin (B) tr(A) for positive semi-definite A, B, GW22 (X, Y) ≥ 4

n X

υi tr(Σi Σξ ) ≥ 4λmin (Σξ )

i=1

= 4λmin (Σξ )

n X

υi tr(Σi )

i=1 n Z X i=1

2 ∥x − mi ∥2 dπi (x) ≥ 4λmin (Σξ )qW n (ξ) ,

X

where the last step follows from (2.2). Taking square roots and the infimum over Y yields Item (i). We show Item (ii), i.e. g = ∥· − ·∥. Fix i ∈ [n] and y ∈ X. Since ∥x − y∥ + ∥x′ − y∥ ≤ 2 diam(X) for x, x′ ∈ X, it holds
2
2
2 ∥x − y∥2 − ∥x′ − y∥2 = ∥x − y∥ − ∥x′ − y∥ ∥x − y∥ + ∥x′ − y∥

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

35

2 ≤ 4 diam(X)2 ∥x − y∥ − ∥x′ − y∥ . Integrating over (x, x′ ) with respect to π̄i ⊗ π̄i and using the second equality in (C.2), we get Varπ̄i (∥· − y∥2 ) ≤ 4 diam(X)2 Varπ̄i (∥· − y∥). Together with (C.3) and the estimate from Item (i), we arrive at Z n X
λmin (Σξ ) W 2 1 2 2 Var ∥· − y∥ dξ(y) ≥ q (ξ) . υ GW2 (X, Y) ≥ π̄ i i 2 4 diam(X) diam(X)2 n X i=1

Taking square roots and the infimum over Y yields Item (ii). C.3. Two-Sided Bounds for the Scalar-Product Gauge. For the proof of Theorem 4.5, we require two auxiliary lemmas. The first one ensures that for any X = (X, ⟨·, ·⟩, ξ) ∈ GM with X ⊂ Rd , the measure ξ has finite second moment. Lemma C.3. Let X ⊂ Rd and ξ ∈ P(X). If ⟨·, ·⟩ ∈ L2 (ξ ⊗ ξ), then ξ ∈ P (2) (X). RR Proof. Set I := Rd ×Rd ⟨x, x′ ⟩2 dξ(x) dξ(x′ ) < ∞. For R > 0, let BR := {x ∈ X : ∥x∥ ≤ R} and Z MR := xx⊤ dξ(x) ∈ Rd×d , BR

which is well-defined and symmetric positive semi-definite, with entries bounded by R2 . Let λ1 , . . . , λd ≥ 0 be the eigenvalues of MR . Then the Cauchy-Schwarz inequality yields X 2 X d d d d X  X 2 2 2 tr(MR ) = λi ≤ λi λ2i = d ∥MR ∥2F . 1 =d i=1

i=1

i=1

i=1

Furthermore, due to tr(xx⊤ ) = ∥x∥2 and ⟨a, b⟩2 = tr(aa⊤ bb⊤ ), we obtain Z 2 2 Z 2 ∥x∥ dξ(x) = tr xx⊤ dξ(x) = tr(MR )2 ≤ d ∥MR ∥2F = d tr(MR MR ) BR Z ZBR ZZ ⊤ ′ ′ ⊤ ′ =d tr(xx x (x ) ) dξ(x) dξ(x ) = d ⟨x, x′ ⟩2 dξ(x) dξ(x′ ) ≤ d I. 2 BR

2 BR

√ R Hence BR ∥x∥2 dξ ≤ d I for every R > 0, and taking the limit R → ∞, monotone convergence yields Z √ ∥x∥2 dξ(x) ≤ d I < ∞, X

as desired. For the remainder of the section, for π ∈ Π(ξ, ∗n ), we define Z n X ai (π)(ai (π))⊤ d (C.4) ai (π) := x dπi (x) ∈ R , i ∈ [n] and M (π) := ∈ Rd×d , πi (X) X i=1 R where we set 0/0 = 0. This makes M (π) well-defined as ∥ai (π)∥2 ≤ πi (X) X ∥x∥2 dπi (x) by Jensen’s inequality, so that ai (π) = 0 whenever πi (X) = 0. In the scalar product case, the quantization problem is equivalent to maximization of the function π 7→ ∥M (π)∥2F , which is the result of the following lemma.

36

F. BEIER, S. ECKSTEIN

Lemma C.4. Let X = (X, ⟨·, ·⟩, ξ) ∈ GM with X ⊂ Rd . Let Mξ := π ∈ Π(ξ, ∗n )
inf QGW G, (πi )ni=1 = ∥Mξ ∥2F − ∥M (π)∥2F .

R

X xx

⊤ dξ.

For every

(n)

G∈Gsym

Proof. Note that Mξ is well-defined by Lemma C.3. Indices i ∈ [n] with πi (X) = 0 contribute to neither side of the claimed identity, so we may assume πi (X) > 0 for all i ∈ [n]. For fixed π ∈ Π(ξ, ∗n ), the minimization of G 7→ QGW (G, (πi )ni=1 ) is solved by Ĝ from Theorem 3.2. Due to the bilinearity of the scalar product it is given by ZZ ⟨ai (π), ai′ (π)⟩ 1 ⟨x, x′ ⟩ dπi (x) dπi′ (x′ ) = Ĝi,i′ = , i, i′ ∈ [n]. πi (X)πi′ (X) X 2 πi (X)πi′ (X) Applying Lemma A.2 with a0 = 0 as usual yields inf (n)

QGW (G, (πi )ni=1 ) =

ZZ X2

G∈Gsym

n X ⟨ai (π), ai′ (π)⟩2

⟨x, x′ ⟩2 dξ(x) dξ(x′ ) −

i,i′ =1

πi (X)πi′ (X)

.

Due to ⟨a, b⟩2 = tr(aa⊤ bb⊤ ) for a, b ∈ Rd , the first term equals ∥Mξ ∥2F and the second term equals tr(M (π)2 ) = ∥M (π)∥2F , as desired. The previous result allows us to derive appropriate bounds in the scalar product case. Before we proceed with this, we recall the Löwner order of symmetric matrices. That is, for two symmetric matrices A, B, we write A ⪰ B if A − B is positive semi-definite. In particular, A is positive semi-definite if and only if A ⪰ 0. We also write A ⪯ B if B − A is positive semi-definite. The bounds rely on the trace/covariance structure of the scalar product GW problem, which has also been exploited e.g. in [23, 69, 78]. Proof of Theorem 4.5. By Lemma C.3, it holds ξ ∈ P (2) (X), so qW n (ξ) and Mξ are welldefined. Before we show the desired bounds, we express qGW (X) via a matrix trace. For n π ∈ Π(ξ, ∗n ) consider ai (π), i ∈ [n], as well as M (π) defined in (C.4). By (3.3) and Lemma C.4,
2 qGW inf inf QGW (G, (πi )ni=1 ) = inf ∥Mξ ∥2F − ∥M (π)∥2F . n (X) = π∈Π(ξ,∗n ) G∈G (n)

π∈Π(ξ,∗n )

sym

Now further define M̃ (π) := Mξ − M (π). As Mξ and M (π) are symmetric, for every π ∈ Π(ξ, ∗n )
∥Mξ ∥2F − ∥M (π)∥2F = tr(Mξ2 ) − tr(M (π)2 ) = tr M̃ (π) (Mξ + M (π)) , hence (C.5)

2 qGW n (X) =

inf π∈Π(ξ,∗n )

tr M̃ (π) (Mξ + M (π)) .

We remark that it holds M̃ (π) =

n Z  X i=1

X

x−

ai (π)  ai (π) ⊤ dπi (x), x− πi (X) πi (X)

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

37

P which follows by expanding the integrand and using ni=1 πi = ξ. We note that M̃ (π) as well as M (π) = Mξ − M̃ (π) are positive semi-definite, hence Mξ ⪯ Mξ + M (π) ⪯ 2Mξ

so that

λmin (Mξ )I ⪯ Mξ + M (π) ⪯ 2λmax (Mξ )I.

Since M̃ (π) ⪰ 0, the map A 7→ tr(M̃ (π) A) is monotone with respect to the Löwner order ⪯. Combining this with (C.5) and the monotonicity of the infimum, we obtain
λmin (Mξ ) inf tr(M̃ (π)) ≤ inf tr M̃ (π) (Mξ + M (π)) π∈Π(ξ,∗n ) π∈Π(ξ,∗n ) | {z } 2 =qGW n (X)

≤ 2λmax (Mξ )

(C.6)

inf

tr(M̃ (π)).

π∈Π(ξ,∗n )

2 Finally, note that inf π∈Π(ξ,∗n ) tr(M̃ (π)) = qW n (ξ) . Indeed, n Z n Z  2  1 Z X X ai (π) 2 tr(M̃ (π)) = x− dπi (x) dπi (x) = x′ dπi (x′ ) x− πi (X) πi (X) X i=1 X i=1 X n Z X = inf ∥x − yi ∥2 dπi (x), n (yi )i=1 ⊂X

i=1

X

where the last line follows by Lemma A.2 and closed convexity of X ⊂ Rd . Taking the infimum over π ∈ Π(ξ, ∗n ) and comparing with (2.2) gives the claim. Plugging this into (C.6) and taking the square root yields the desired bounds. The smallest eigenvalue in the lower bound of Theorem 4.5 may render it unnecessarily weak, particularly when it is zero. A sharper estimate can be achieved by applying the following result. Lemma C.5 (Ruhe trace inequality including subspaces). Let A, B ∈ Rn×n be symmetric and positive semi-definite with range(A) ⊆ range(B) =: S and s := dim(S). Let λ1 (B) ≥ . . . ≥ λs (B) > 0 be the nonzero eigenvalues of B and λ1 (A) ≥ . . . ≥ λs (A) ≥ 0 the s largest eigenvalues of A (which can include zeros). Then tr(AB) ≥

s X

λi (A)λs+1−i (B).

i=1

Proof. By symmetry of A and range(A) ⊆ S, we have S ⊥ ⊆ ker(A), so S is invariant under A, and analogously under B. Hence tr(AB) = tr(A|S B|S ). Applying the standard Ruhe trace inequality (see [40, 9.H.1.h.]) yields the claim. The proof of Theorem 4.5 implicitly relies on the estimate tr(M̃ Mξ ) ≥ λmin (Mξ ) tr(M̃ ) √ which amounts to scaling ξ by λmin along every direction. Since 0 ⪯ M̃ ⪯ Mξ implies range(M̃ ) ⊆ range(Mξ ), applying Lemma C.5 instead gives the sharper estimate tr(M̃ Mξ ) ≥

s X

λi (M̃ ) λs+1−i (Mξ ),

s := rank(Mξ ).

i=1

This scales the directions of M̃ by the varying nonzero eigenvalues of Mξ . In particular, λmin in Theorem 4.5 may be replaced by the smallest nonzero eigenvalue of Mξ .

38

F. BEIER, S. ECKSTEIN

C.4. Monge Maps for the Scalar-Product Gauge. In the following, we prove Theorem 4.7. Due to (3.3) and Lemma C.4 the quantization problem qGW n (X) reduces to the 2 problem supπ∈Π(ξ,∗n ) ∥M (π)∥F . We use this reduction in the following to establish the existence of Monge maps from X to solutions of the GW quantizers. By Theorem 3.4, this has direct implications on the structure of the solutions themselves and gives a reformulation of the problem that is akin to the partitioning problem (2.4). Before proving Theorem 4.7, we establish some fundamental properties of the reduced objective in the following lemma. The proofs of both the lemma and the theorem rely on the notation ai (π), i ∈ [n], and M (π), π ∈ Π(ξ, ∗n ), introduced in (C.4). Lemma C.6. The map π 7→ ∥M (π)∥2F is convex and weakly continuous on Π(ξ, ∗n ). Proof. Set F : Π(ξ, ∗n ) → [0, ∞), π 7→ ∥M (π)∥2F . For convexity, let π (0) , π (1) ∈ Π(ξ, ∗n ), t ∈ [0, 1] and π (t) := (1 − t)π (0) + tπ (1) ∈ Π(ξ, ∗n ). Note that (ai (π), πi (X)) depends linearly on π. Furthermore, the joint convexity of the matrixfractional function (a, p) 7→ aa⊤ /p in the Löwner order (see e.g. [14, Sec. 3.1.7]) gives, for each (0) (1) i ∈ [n] with πi (X), πi (X) > 0, ai (π (t) )(ai (π (t) ))⊤ (t)

πi (X)

⪯ (1 − t)

ai (π (0) )(ai (π (0) ))⊤ (0)

+t

ai (π (1) )(ai (π (1) ))⊤ (1)

πi (X)

.

πi (X) (0)

(1)

For the remaining indices the inequality holds directly: If πi (X) = πi (X) = 0, then ai (π (0) ) = ai (π (1) ) = ai (π (t) ) = 0 and all three terms vanish. If exactly one vanishes, say (0) (t) (1) πi (X) = 0 (so ai (π (0) ) = 0), then by linearity ai (π (t) ) = t ai (π (1) ) and πi (X) = t πi (X). (1) Thus, for t ∈ (0, 1) both sides equal t ai (π (1) )(ai (π (1) ))⊤ /πi (X) (the cases t ∈ {0, 1} being trivial). Hence the inequality holds for every i ∈ [n], and summing over i ∈ [n] yields 0 ⪯ M (π (t) ) ⪯ (1 − t) M (π (0) ) + t M (π (1) ). Moreover, ∥·∥F is monotone with respect to ⪯ for positive semi-definite matrices, i.e. for 0 ⪯ A ⪯ B, ∥A∥2F = tr(A2 ) ≤ tr(AB) ≤ tr(B 2 ) = ∥B∥2F , as the trace of a product of two positive semi-definite matrices is non-negative. Combining both with the triangle inequality and the convexity of s 7→ s2 , we obtain 

2  2 F (π (t) ) = tr M (π (t) )2 ≤ tr (1 − t)M (π (0) ) + tM (π (1) ) = (1 − t)M (π (0) ) + tM (π (1) ) F ≤ (1 − t)∥M (π (0) )∥F + t∥M (π (1) )∥F

2

≤ (1 − t)F (π (0) ) + tF (π (1) ),

that is, F is convex. For continuity, let π (k) ⇀ π in Π(ξ, ∗n ), we show M (π (k) ) → M (π), k → ∞. Fix a continuous cutoff χR : X → [0, 1] with 1{∥x∥≤R} ≤ χR (x) ≤ 1{∥x∥≤R+1} , x ∈ X. Each coordinate of

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING (k)

39

R (k) → X x χR dπi as X x χR dπi (k) k → ∞. Moreover, since 1 − χR is supported on {∥x∥ > R} and by πi ≤ ξ, we obtain x 7→ x χR (x) is bounded and continuous, so πi Z X

Z

(k) x (1 − χR ) dπi

≤ {∥x∥>R}

⇀ πi gives

(k) ∥x∥ dπi

R

Z ≤

∥x∥ dξ → 0, {∥x∥>R} (k)

uniformly in k, and the same bound holds with πi in place of πi . Hence, for every R > 0, Z Z Z Z (k) (k) (k) ∥x∥ dπi . ∥x∥ dπi + x χR dπi + x χR dπi − ai (π ) − ai (π) ≤ {∥x∥>R}

{∥x∥>R}

X

X

Letting first k → ∞, the first term vanishes while the remaining terms are each bounded by R (k) ) → a (π). If π (X) > 0, the i-th summand i i {∥x∥>R} ∥x∥ dξ. Letting then R → ∞ gives ai (π of M (π (k) ) converges to that of M (π). If πi (X) = 0, then for any R > 0, Jensen’s inequality yields Z Z ai (π (k) )(ai (π (k) ))⊤ ∥ai (π (k) )∥2 (k) 2 2 (k) = ≤ ∥x∥ dπi ≤ R πi (X) + ∥x∥2 dξ. (k) (k) F | {z } X {∥x∥>R} π (X) π (X) i

i

→0,k→∞

Hence, lim sup

ai (π (k) )(ai (π (k) ))⊤ (k) πi (X)

k→∞

Z F

≤

R→∞

∥x∥2 dξ −−−−→ 0.

{∥x∥>R}

Thus, M (π (k) ) → M (π) and F is weakly continuous on Π(ξ, ∗n ). Proof of Theorem 4.7. Using (3.3) together with Lemma C.4, the quantization problem (3.2) becomes equivalent to maximizing F (π) := ∥M (π)∥2F over Π(ξ, ∗n ). By Lemma C.6, F is convex and weakly continuous on Π(ξ, ∗n ). Moreover, Π(ξ, ∗n ) is convex and weakly compact. Thus, the Bauer maximum principle, see e.g. [2, Thm. 7.69], ensures that F attains its maximum over Π(ξ, ∗n ) at an extreme point of Π(ξ, ∗n ). Now, we show that every extreme point of Π(ξ, ∗n ) is induced by a map. Let π ∈ Π(ξ, ∗n ) R be extreme and let π = X δx ⊗ κx dξ(x) be its disintegration, that is Z Z Z f (x, i) dπ(x, i) = f (x, i) dκx (i) dξ(x) for all bounded measurable f. X×[n]

X

[n]

To achieve a contradiction, we assume that B := {x ∈ X : maxi∈[n] κx ({i}) < 1} satisfies ξ(B) > 0. For x ∈ B, let ix < jx denote the two smallest indices with κx ({ix }) > 0, κx ({jx }) > 0 and define the measurable family of signed measures 
σx := min κx ({ix }), κx ({jx }) · δix − δjx 1B (x). R Clearly σx ([n]) = 0 for all x ∈ X, hence π ± := X δx ⊗ (κx ± σx ) dξ(x) define two distinct elements of Π(ξ, ∗n ). Furthermore, π = 21 (π + + π − ), contradicting extremality. Hence ξ(B) =

40

F. BEIER, S. ECKSTEIN

0 and thus κx is a Dirac measure for ξ-a.e. x ∈ X, i.e., π = (Id, T )# ξ for a measurable map T : X → [n]. By the first part of the proof, there exists a measurable T : X → [n] such that π ∗ := (Id, T )# ξ maximizes F over Π(ξ, ∗n ) and thus also minimizes (C.7)

inf (n)

QGW (G, (πi )ni=1 ) = ∥Mξ ∥2F − ∥M (π)∥2F

G∈Gsym 2 ∗ n over Π(ξ, ∗n ), cf. Lemma C.4. By (3.3), this implies qGW n (X) = inf G∈G (n) QGW (G, (πi )i=1 ). sym

We set Vi := T −1 ({i}). For indices i ∈ [n] with ξ(Vi ) = 0, it holds πi∗ = 0, so they contribute to neither side of (C.7), and the corresponding entries of G do not affect QGW (G, (πi∗ )ni=1 ). Hence we may apply Theorem 3.2 to the indices with ξ(Vi ) > 0. Thus, the latter is minimized by the block-mean gauge, which by bilinearity of the scalar product equals     Z Z ai (π ∗ ) ai′ (π ∗ ) 1 1 Ĝi,i′ = = , x dξ(x), x dξ(x) = ⟨mξ (Vi ), mξ (Vi′ )⟩, πi∗ (X) πi∗′ (X) ξ(Vi ) Vi ξ(Vi′ ) Vi′ where we used πi∗ = ξ|Vi . Now set υ̂ = As π ∗ ∈ Π(ξ, υ̂),

Pn

i=1 ξ(Vi )δi giving rise to the gm-space Ŷ = ([n], Ĝ, υ̂).

2 2 ∗ GW 2 qGW n (X) ≤ GW2 (X, Ŷ) ≤ QGW (Ĝ, (πi )) = qn (X) ,

so Ŷ is a solution of the quantization problem and π ∗ = (Id, T )# ξ an optimal coupling. The crux of the argument is that every extreme point of Π(ξ, ∗n ) is induced by a map. We establish this directly, by a perturbation argument tailored to our setting. Alternatively, it can be R derived from the theory of Young measures. More precisely, under the disintegration π = X δx ⊗κx dξ(x), the set Π(ξ, ∗n ) is identified with the transition kernels x 7→ κx ∈ P([n]), whose extreme points are precisely the Dirac-valued ones (see e.g. [17]). C.5. Closed-Form for the Scalar-Product Case in 1D. Proof of Theorem 4.8. As ∥u∥ = 1, we obtain ⊤ ′ ⊤ ′ ′ T (x) T (x′ ) = u⊤ x u⊤ x′ = x⊤ uu | {z x} = x x = ⟨x, x ⟩,

x, x′ ∈ X ⊂ span(u).

=x′

˜ For a partition Thus, setting ξ˜ := T# ξ we see that X is homomorphic to X̃ = (R, (s, t) 7→ st, ξ). n (Vi )i=1 of R write n X
˜ i ) m (Vi )2 . ξ(V S (Vi )ni=1 := ξ̃ i=1 GW Since X and X̃ are homomorphic, it holds qGW n (X) = qn (X̃). Using (4.3), we obtain 2 GW 2 qGW n (X) = qn (X̃) =

infn

(Vi )i=1

n Z X i,i′ =1 Vi ×Vi′

st − mξ̃ (Vi )mξ̃ (Vi′ )

2

˜ dξ(t) ˜ dξ(s)

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

= infn

(Vi )i=1

= infn

(Vi )i=1

= infn

(Vi )i=1

n Z X Vi ×Vi′

i,i′ =1

Z R

41

 ˜ i )ξ(V ˜ i′ ) ˜ dξ(t) ˜ − m (Vi )2 m (Vi′ )2 ξ(V (st)2 dξ(s) ξ̃ ξ̃

n 2  2  X ˜ ˜ i ) m (Vi )2 t dξ − ξ(V ξ̃ 2

i=1

˜ 2 − S((Vi ))2 = M2 (ξ) ˜ 2 − sup S((Vi ))2 , M2 (ξ) (Vi )

R ˜ i ) m (Vi ). On the other hand, a similar where equality in the second line uses Vi s dξ˜ = ξ(V ξ̃ derivation yields Z X ˜ = M2 (ξ) ˜ − S((Vi )), (t − mξ̃ (Vi ))2 dξ(t) i

Vi

n ˜2 ˜ so that (2.4) gives qW n (ξ) = M2 (ξ)−sup(Vi ) S((Vi )i=1 ). Both problems amount to maximizing
2 S, and since S ≥ 0 it holds sup(Vi ) S((Vi ))2 = sup(Vi ) S((Vi )) . Therefore

2 2 W ˜2 2 W ˜2 ˜2 ˜2 ˜ ˜2 ˜ qGW = qW n (X) = M2 (ξ) −sup S((Vi )) = M2 (ξ) − M2 (ξ)−qn (ξ) n (ξ) 2M2 (ξ)−qn (ξ) , (Vi )

as desired. Appendix D. Proofs of Section 5. D.1. Directional Derivative of R. Proof of Proposition 5.1. Let η ∈ Π(ξ, ∗n ). Write πt := (1 − t)π + tη = π + t(η − π). Since ηi ∈ M+ (X) and πi (X) > 0 by assumption, it holds (πt )i (X) = (1 − t)πi (X) + t ηi (X) > 0, so Ĝ(πt ) is well-defined. Let ZZ ai,i′ (t) := g d(πt )i d(πt )i′ ,

t ∈ [0, 1)

bi (t) := (πt )i (X) = (1 − t) πi (X) + t ηi (X),

X2 a

′ (t)

i,i so that Ĝ(πt )i,i′ = bi (t)b . Since bi (t) = (1 − t) πi (X) + t ηi (X) ≥ (1 − t) πi (X) > 0 for ′ (t) i

t ∈ [0, 1), each Ĝ(πt )i,i′ is a quotient of polynomials in t with non-vanishing denominator, hence differentiable on [0, 1). Similarly, as t 7→ πt is affine, X ZZ
2 t 7→ QGW (G, πt ) = g − Gi,i′ d(πt )i d(πt )i′ i,i′ ∈[n]

X2

is a quadratic polynomial. In addition, QGW is a quadratic polynomial in its first argument. Therefore, as a composition of differentiable functions, t 7→ R(πt ) = QGW (Ĝ(πt ), πt ) exists and equals is differentiable on [0, 1). Hence, the one-sided limit limt↓0 R(πt )−R(π) t d dt t=0 R(πt ), so it suffices to compute the latter. As Ĝ(πt ) is the unconstrained minimizer

42

F. BEIER, S. ECKSTEIN

of the strictly convex quadratic G 7→ QGW (G, πt ) (Theorem 3.2), its first-order optimality gives ∂G QGW (Ĝ(πt ), πt ) = 0, so the chain rule gives d  d dπt Ĝ(πt ) + ∂π QGW (Ĝ(πt ), πt )[ ]. R(πt ) = ∂G QGW (Ĝ(πt ), πt ) dt {z } | dt dt |{z} =0

=(η−π)

For the remaining derivative, we see that QGW (G, γ + s α) − QGW (G, γ) s X ZZ
(g − Gi,i′ )2 dαi dγi′ + dγi dαi′ + s dαi dαi′ = lim

∂π QGW (G, γ)[α] = lim

s→0

s→0

=

X2

i,i′ ∈[n]

X ZZ

(g − Gi,i′ )2 dαi dγi′ + dγi dαi′

X2

i,i′ ∈[n]

X ZZ

=2

i,i′ ∈[n]

(g − Gi,i′ )2 dαi dγi′ , X2

(n) G ∈ Gsym , γ ∈ Π(ξ, ∗n ),

where the last line follows by symmetry of g and G. In summary, we have X ZZ
2 d R(πt ) = ∂π QGW (Ĝ(πt ), πt )[η − π] = 2 g − Ĝ(πt )i,i′ d(πt )i′ d(ηi − πi ). dt X2 ′ i,i ∈[n]

The formula holds for all t ∈ [0, 1), in particular for t = 0, we obtain n Z X ZZ X
2 d (π) R(πt ) = 2 g − Ĝ(π)i,i′ dπi′ d(ηi − πi ) = 2 ĉi d(ηi − πi ), dt t=0 2 X X ′ i=1

i,i ∈[n]

as desired. D.2. Convergence to Stationary Quantizers. Proof of Theorem 5.3. Let (π (k) )k be generated by Algorithm 1 with the line search of (k) Remark 5.2, i.e. π (k+1) = (1 − tk )π (k) + tk η (k) with ηi = ξ|V (k) for a Voronoi partition i
(k) (π (k) ) n (Vi )ni=1 with respect to (ĉi )i=1 , and tk ∈ argmins∈[0,1] QGW Ĝ(π (k) ), (1 − s)π (k) + sη (k) . Analogously to Remark 5.2, let ak and bk be given by n ZZ X
2 (k) (k) (k) (k) ak := g − Ĝ(π (k) )i,i′ d(ηi − πi ) d(ηi′ − πi′ ), X2 i,i′ =1 n Z X

bk := −2

i=1

X

(π (k) )

ĉi

(k)

d(ηi

(k)

− πi ).

P R (π(k) ) Note that bk is non-negative since, by Lemma A.1, η (k) minimizes η 7→ i X ĉi dηi over Π(ξ, ∗n ). Before showing stationarity, we show that bk → 0, as k → ∞. It holds
QGW Ĝ(π (k) ), (1 − s)π (k) + sη (k) = R(π (k) ) − s bk + s2 ak , s ∈ [0, 1].

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

43

By construction, tk minimizes this quadratic over [0, 1], so that, with π (k+1) = (1 − tk )π (k) + tk η (k) ,

QGW Ĝ(π (k) ), π (k+1) = min R(π (k) ) − s bk + s2 ak . s∈[0,1]

Since R(π) = inf G QGW (G, π) ≤ QGW (Ĝ(π (k) ), π) for every π ∈ Π(ξ, ∗n ), in particular R(π (k+1) ) ≤ QGW (Ĝ(π (k) ), π (k+1) ), and therefore

R(π (k) ) − R(π (k+1) ) ≥ R(π (k) ) − min R(π (k) ) − s bk + s2 ak = max s bk − s2 ak . s∈[0,1]

s∈[0,1]

Each Ĝ(π (k) )i,i′ is an average of g and since g is bounded, we have |Ĝ(π (k) )i,i′ | ≤ ∥g∥∞ and (k) (k) thus (g − Ĝ(π (k) )i,i′ )2 ≤ 4∥g∥2∞ . As each signed measure ηi − πi has total variation at most 2, it follows that ak is bounded above by some A ∈ (0, ∞) uniformly in k. Using this bound, we get o n
b2k , R(π (k) ) − R(π (k+1) ) ≥ max s bk − s2 A ≥ min b2k , 4A s∈[0,1]

where the last inequality holds since the concave function s 7→ s bk − s2 A is maximized on b2k bk bk }, yielding 4A (if s∗ = 2A ) or bk − A ≥ b2k (if s∗ = 1, i.e., bk ≥ 2A). [0, 1] at s∗ = min{1, 2A Summing the per-step bound, for every K ∈ N, we obtain R(π (0) ) ≥ R(π (0) ) − R(π (K) ) =

K−1 X

n o X
K−1 b2k R(π (k) ) − R(π (k+1) ) ≥ min b2k , 4A .

k=0

k=0

 b b2 P (0) k k Letting K → ∞ yields ∞ k=0 min 2 , 4A ≤ R(π ) < ∞, so bk → 0 as k → ∞. Finally, we show the desired stationarity. Let π ∗ be a cluster point with πi∗ (X) > 0 for all i. By compactness of Π(ξ, ∗n ), pass to subsequences π (kj ) ⇀ π ∗ and η (kj ) ⇀ η̄ ∈ Π(ξ, ∗n ). We first record the convergence Z Z (kj ) j→∞ (π (kj ) ) (π ∗ ) (D.1) ĉi dσi −−−→ ĉi dσi , i ∈ [n], X

X

(k )

valid for any (σ (kj ) )j ⊂ Π(ξ, ∗n ) with σ (kj ) ⇀ σ. Indeed, firstly since πi j (X) → πi∗ (X) > 0 and, by Lemma B.1 with h = g, the numerators of Ĝ(π (kj ) )i,i′ converge, we have Ĝ(π (kj ) )i,i′ → Ĝ(π ∗ )i,i′ for all i, i′ ∈ [n]. Expanding the square in Z X

(k ) (π (kj ) ) ĉi dσi j =

n ZZ X i′ =1

g − Ĝ(π (kj ) )i,i′ X2

2

(k )

(k )

dσi j dπi′ j

and applying Lemma B.1 to h ∈ {g 2 , g, 1} ⊂ L1 (ξ ⊗ ξ) (using that g is bounded) yields (D.1). P R (π∗ ) Applying (D.1) with σ (kj ) = η (kj ) and with σ (kj ) = π (kj ) gives bkj → −2 i X ĉi d(η̄i − πi∗ ) as j → ∞. By the first part, (bkj )j∈N also converges to 0, so that X Z (π∗ ) (D.2) 2 ĉi d(η̄i − πi∗ ) = 0. i

X

44

F. BEIER, S. ECKSTEIN

We conclude the proof by showing X Z (π∗ ) X Z (π∗ ) ĉi dη̄i ĉi dηi ≥ 2 (D.3) 2 X

i

for all η ∈ Π(ξ, ∗n ).

X

i

P R (π(k) ) By construction, Lemma A.1 ensures that η (k) minimizes η 7→ 2 i X ĉi dηi over Π(ξ, ∗n ), so for every η ∈ Π(ξ, ∗n ), X Z (π(k) ) (k) X Z (π(k) ) ĉi dηi . ĉi dηi ≥ 2 2 X

i

X

i

We evaluate this inequality along the subsequence (kj )j and pass to the limit via (D.1). P R (π∗ ) Firstly, with the constant sequence σ (kj ) ≡ η the left-hand side converges to 2 i X ĉi dηi . P R (π∗ ) Secondly, with σ (kj ) = η (kj ) the right-hand side converges to 2 i X ĉi dη̄i and we obtain (D.3). Combining this inequality with (D.2), for every η ∈ Π(ξ, ∗n ), X Z (π∗ ) X Z (π∗ ) ∗ ĉi d(η̄i − πi∗ ) = 0, ĉi d(ηi − πi ) ≥ 2 2 i

X

X

i

so π ∗ is stationary, as desired. D.3. Unit-Step Descent for the Scalar-Product Gauge. (k)

Proof of Corollary 5.4. From Algorithm 1, let (π (k) )k∈N and ((Vi )ni=1 )k∈N be the sequence of measures and Voronoi partitions, respectively. By Lemma C.4 and Lemma C.6, R is concave on Π(ξ, ∗n ), that is (1 − t)R(π) + tR(γ) ≤ R((1 − t)π + tγ) which can be rearranged as R(γ) − R(π) ≤

R((1 − t)π + tγ) − R(π) , t

π, γ ∈ Π(ξ, ∗n ), t ∈ (0, 1].

Let k ∈ N and take π = π (k) as well as γ = π (k+1) in the above inequality. Taking the limit t ↓ 0 and applying Proposition 5.1 yields R(π

(k+1)

) − R(π

(k)

)≤2

n Z X i=1

(k+1)

Furthermore, by construction πi

In summary R(π (k+1) ) ≤ R(π (k) ).

(k+1)

d(πi

X

(π (k) )

ĉi

(k+1)

d(πi

(k)

− πi ) ≤ 0.

(k)

− πi ).

= ξ|V (k) , so Lemma A.1 yields i

n Z X i=1

X

(π (k) )

ĉi

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

45

Assume now that X is finite. Then Π(ξ, ∗n ) is a polytope whose vertices are the hard (k+1) assignments πi = ξ|Vi , and the unit-step update πi = ξ|V (k) produces such a vertex at i

(π (k) )

(k)

(π (k) )

is not concentrated on {x ∈ X : ĉi (x) = minj ĉj (x)}, i ∈ [n], then P R (π(k) ) (k+1) (k) by Lemma A.1 the inequality i ĉi d(πi − πi ) ≤ 0 is strict, so that R(π (k+1) ) < (k) R(π ). Since R then strictly decreases over the finitely many vertex values, this can occur only finitely often; hence after finitely many steps π (k) is concentrated on the desired sets. every step. If πi

Appendix E. Non-Uniformity of the Lower Quantization Bound. Although Corollaries 4.4 and 4.6 pin down the rate, the constants c to bound qGW by n W c · qn from below genuinely depend on the measure and cannot be chosen uniformly over GM. A single family witnesses this simultaneously for any of the three gauges g ∈ {∥ · − · ∥, ∥ · − · ∥2 , ⟨·, ·⟩}. More precisely, let X := X(n) := (R2 , g, ξ (n) ), where n−1

ξ (n) :=

X  1  δqj , δp+ + δp− + n+1

p± := (0, ±1), qj := (j, 0).

j=1

We consider a quantization p with n points. Since each qj has vanishing second coordinate, ⟨p± , qj ⟩ = 0 and ∥p± −qj ∥ = j 2 + 1 are independent of the sign, so p+ and p− carry identical gauge to every qj , for each such g. Consider the n-point quantizer obtained by merging p+ , p− into a single (abstract) point z. Then, we set Y = (Y, h, υ), where Y := {z, q1 , . . . , qn−1 },  X 1  g(ps , ps′ ) y = y ′ = z,    4  s,s′ ∈{+,−} h(y, y ′ ) = g(p+ , qj ) {y, y ′ } = {z, qj }, j ∈ [n − 1],      g(qj , qk ) y = qj , y ′ = qk , j, k ∈ [n − 1]. 1 2 δz + n+1 and υ = n+1

Pn−1

j=1 δqj . Under the canonical coupling,

n−1  X
1  δ(p+ , z) + δ(p− , z) + δ(qj , qj ) ∈ Π ξ (n) , υ , π= n+1 j=1

all pairs involving some qj contribute zero distortion. Consequently, the only nonzero distortion under π stems from the four pairs in {p+ , p− }2 , each of which is mapped to (z, z). Hence ZZ
2 2 g(x, x′ ) − h(y, y ′ ) dπ(x, y) dπ(x′ , y ′ ) GW2 (X, Y) ≤ =

1 (n + 1)2

X

2 g(ps , ps′ ) − h(z, z) .

s,s′ ∈{+,−}

|

{z

=:Cg

}

The sum in Cg depends only on the fixed points p± = (0, ±1) and not on n, so it is a constant Cg < ∞ for each of the three gauges. Moreover Cg > 0. Indeed, since g(p+ , p− ) ̸= g(p+ , p+ ),

46

F. BEIER, S. ECKSTEIN

the four values g(ps , ps′ ) are not all equal, and hence their deviation from the mean h(z, z) does not vanish. Since Y ∈ GMn , this yields 2 2 qGW n (X) ≤ GW2 (X, Y) ≤

Cg . (n + 1)2

(n) ) is taken with respect to the Euclidean metric For the Wasserstein value, recall that qW n (ξ in all three cases and (cf. (2.3)) Z (n) 2 qW (ξ ) = inf d2 (x, Y ′ ) dξ (n) (x), d(x, Y ′ ) := min′ ∥x − y∥. n |Y ′ |≤n

y∈Y

Fix any Y ′ ⊂ R2 with |Y ′ | ≤ n. Since |supp(ξ (n) )| = n + 1 and |Y ′ | ≤ n, the nearest-point assignment that maps x ∈ supp(ξ (n) ) to a minimizer of Y ′ → [0, ∞), y 7→ ∥x − y∥, is not injective. Thus, we can fix two distinct atoms x0 ̸= x′0 ∈ supp(ξ (n) ) with a common nearest point y0 ∈ Y ′ . As the atoms of ξ (n) are pairwise at Euclidean distance at least 1, it holds 1 ≤ ∥x0 − x′0 ∥ ≤ ∥x0 − y0 ∥ + ∥y0 − x′0 ∥ ≤ 2 max{d(x0 , Y ′ ), d(x′0 , Y ′ )}. Hence, at least one of them, say x0 , satisfies d(x0 , Y ′ ) ≥ 21 . Keeping only this atom’s contribution to the integral, Z 1 1 · . d2 (x, Y ′ ) dξ (n) (x) ≥ ξ (n) ({x0 }) d2 (x0 , Y ′ ) ≥ n+1 4 1 (n) )2 ≥ Since Y ′ was arbitrary, qW n (ξ 4(n+1) . Combining the two bounds,

p p 2 Cg n→∞ Cg /(n + 1)2 qGW n (X) =√ ≤p −−−→ 0. (n) ) qW n+1 1/(4(n + 1)) n (ξ Hence no constant c > 0 with qGW ≥ c qW n n can hold uniformly over GM, for any of the three gauges. AI Tool Disclosure. We used AI, primarily Claude Opus 5 and Fable 5, as follows: First, to identify many inaccuracies in earlier versions of this paper. Second, to generate most of the code for the experiments; we only reviewed the correctness via structural tests and by spot-checking parts of the code. Third, many ideas for the proofs of the lower bounds in Section 4 were suggested by AI. Fourth, many relevant references were located by AI. Finally, AI generated first suggestions for the proofs of some auxiliary lemmas. Aside from the aforementioned suggestions, the text in this paper was human-written. REFERENCES [1] A. K. Akash, S. Li, and N. Garcı́a Trillos, Wasserstein barycenter-based model fusion and linear mode connectivity of neural networks, 2022, https://arxiv.org/abs/2210.06671. arXiv:2210.06671. [2] C. D. Aliprantis and K. C. Border, Infinite Dimensional Analysis: A Hitchhiker’s Guide, Springer, Berlin, 3rd ed., 2006, https://doi.org/10.1007/3-540-29587-9.

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

47

[3] D. Aloise, A. Deshpande, P. Hansen, and P. Popat, NP-hardness of Euclidean sum-of-squares clustering, Mach. Learn., 75 (2009), pp. 245–248, https://doi.org/10.1007/s10994-009-5103-0. [4] D. Alvarez-Melis and T. Jaakkola, Gromov–Wasserstein alignment of word embedding spaces, in Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing, Brussels, Belgium, 2018, Association for Computational Linguistics, pp. 1881–1890, https://doi.org/10.18653/ v1/D18-1214. [5] D. Arthur and S. Vassilvitskii, k-means++: The advantages of careful seeding, in Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA 2007), Philadelphia, PA, 2007, SIAM, pp. 1027–1035. [6] A. D. Aydın, Asymptotics of the quantization problem on metric measure spaces, Math. Ann., 394 (2026), p. 48, https://doi.org/10.1007/s00208-026-03376-x. [7] F. Beier, Gromov–Wasserstein transfer operators, in Scale Space and Variational Methods in Computer Vision (SSVM), vol. 14009 of Lecture Notes in Computer Science, Springer, 2023, pp. 614–626, https://doi.org/10.1007/978-3-031-31975-4 47. [8] F. Beier and R. Beinert, Tangential fixpoint iterations for Gromov–Wasserstein barycenters, SIAM J. Imaging Sci., 18 (2025), pp. 1058–1100, https://doi.org/10.1137/24M1654804. [9] F. Beier, R. Beinert, and G. Steidl, On a linear Gromov–Wasserstein distance, IEEE Trans. Image Process., 31 (2022), pp. 7292–7305, https://doi.org/10.1109/TIP.2022.3221286. [10] F. Beier, R. Beinert, and G. Steidl, Multi-marginal Gromov–Wasserstein transport and barycentres, Inf. Inference, 12 (2023), pp. 2753–2781, https://doi.org/10.1093/imaiai/iaad041. [11] R. Beinert, C. Heiss, and G. Steidl, On assignment problems related to Gromov–Wasserstein distances on the real line, SIAM J. Imaging Sci., 16 (2023), pp. 1028–1032, https://doi.org/10.1137/ 22M1497808. [12] V. I. Bogachev, Weak Convergence of Measures, vol. 234 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, 2018, https://doi.org/10.1090/surv/234. [13] F. Bogo, J. Romero, M. Loper, and M. J. Black, FAUST: Dataset and evaluation for 3D mesh registration, in 2014 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), IEEE, 2014, pp. 3794–3801, https://doi.org/10.1109/CVPR.2014.491. [14] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, Cambridge, 2004, https://doi.org/10.1017/CBO9780511804441. [15] G. Bravo-Hermsdorff and L. M. Gunderson, A unifying framework for spectrum-preserving graph sparsification and coarsening, in Advances in Neural Information Processing Systems, vol. 32, Curran Associates, Inc., 2019, pp. 7734–7745. [16] G. D. Cañas and L. Rosasco, Learning probability measures with respect to optimal transport metrics, in Advances in Neural Information Processing Systems, vol. 25, Curran Associates, Inc., 2012, pp. 2501– 2509. [17] C. Castaing and M. Valadier, Convex Analysis and Measurable Multifunctions, vol. 580 of Lecture Notes in Mathematics, Springer, Berlin, Heidelberg, 1977, https://doi.org/10.1007/BFb0087685. [18] Y. Chen, R. Yao, Y. Yang, and J. Chen, A Gromov–Wasserstein geometric view of spectrumpreserving graph coarsening, in International Conference on Machine Learning (ICML), PMLR, 2023, pp. 5257–5281, https://proceedings.mlr.press/v202/chen23ak.html. [19] S. Chowdhury, D. Miller, and T. Needham, Quantized Gromov–Wasserstein, in Machine Learning and Knowledge Discovery in Databases (ECML PKDD 2021), vol. 12977 of Lecture Notes in Computer Science, Springer, 2021, pp. 811–827, https://doi.org/10.1007/978-3-030-86523-8 49. [20] S. Chowdhury and T. Needham, Gromov–Wasserstein averaging in a Riemannian framework, in 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), IEEE, 2020, pp. 3676–3684, https://doi.org/10.1109/CVPRW50498.2020.00429. [21] S. Chowdhury and T. Needham, Generalized spectral clustering via Gromov–Wasserstein learning, in Proceedings of the 24th International Conference on Artificial Intelligence and Statistics, vol. 130 of Proceedings of Machine Learning Research, PMLR, 2021, pp. 712–720. [22] M. Cuturi and A. Doucet, Fast computation of Wasserstein barycenters, in Proceedings of the 31st International Conference on Machine Learning, E. P. Xing and T. Jebara, eds., vol. 32 of Proceedings of Machine Learning Research, PMLR, 2014, pp. 685–693, https://proceedings.mlr.press/v32/ cuturi14.html.

48

F. BEIER, S. ECKSTEIN

[23] J. Delon, A. Desolneux, and A. Salmona, Gromov–Wasserstein distances between Gaussian distributions, J. Appl. Probab., 59 (2022), pp. 1178–1198, https://doi.org/10.1017/jpr.2022.16. [24] P. Demetci, R. Santorella, B. Sandstede, W. S. Noble, and R. Singh, SCOT: Single-cell multiomics alignment with optimal transport, J. Comput. Biol., 29 (2022), pp. 3–18, https://doi.org/10. 1089/cmb.2021.0446. [25] T. Ding, D. Xiang, P. Rivas, and L. Dong, Neural pruning for 3D scene reconstruction: Efficient NeRF acceleration, in AI Revolution: Research, Ethics and Society, vol. 2723 of Communications in Computer and Information Science, Cham, 2026, Springer, pp. 81–93, https://doi.org/10.1007/ 978-3-032-13056-3 8. [26] T. Dumont, T. Lacombe, and F.-X. Vialard, On the existence of Monge maps for the Gromov– Wasserstein problem, Found. Comput. Math., 25 (2025), pp. 463–510, https://doi.org/10.1007/ s10208-024-09643-0. [27] M. Ehler, M. Gräf, S. Neumayer, and G. Steidl, Curve based approximation of measures on manifolds by discrepancy minimization, Found. Comput. Math., 21 (2021), pp. 1595–1642, https: //doi.org/10.1007/s10208-021-09491-2. [28] R. Flamary, N. Courty, A. Gramfort, M. Z. Alaya, A. Boisbunon, S. Chambon, L. Chapel, A. Corenflos, K. Fatras, N. Fournier, L. Gautheron, N. T. H. Gayraud, H. Janati, A. Rakotomamonjy, I. Redko, A. Rolet, A. Schutz, V. Seguy, D. J. Sutherland, R. Tavenard, A. Tong, and T. Vayer, POT: Python optimal transport, J. Mach. Learn. Res., 22 (2021), pp. 1–8, https://jmlr.org/papers/v22/20-451.html. [29] M. Gräf, D. Potts, and G. Steidl, Quadrature errors, discrepancies, and their relations to halftoning on the torus and the sphere, SIAM J. Sci. Comput., 34 (2012), pp. A2760–A2791, https://doi.org/10. 1137/100814731. [30] S. Graf and H. Luschgy, Foundations of Quantization for Probability Distributions, vol. 1730 of Lecture Notes in Mathematics, Springer, Berlin, 2000, https://doi.org/10.1007/BFb0103945. [31] R. M. Gray and D. L. Neuhoff, Quantization, IEEE Trans. Inf. Theory, 44 (1998), pp. 2325–2383, https://doi.org/10.1109/18.720541. [32] T. Hoefler, D. Alistarh, T. Ben-Nun, N. Dryden, and A. Peste, Sparsity in deep learning: Pruning and growth for efficient inference and training in neural networks, J. Mach. Learn. Res., 22 (2021), pp. 1–124, https://jmlr.org/papers/v22/21-0366.html. [33] G. Houry, J. Feydy, and F.-X. Vialard, Gromov–Wasserstein at scale, beyond squared norms, in International Conference on Machine Learning (ICML), 2026, https://openreview.net/forum?id= P9LofJB8fs, https://arxiv.org/abs/2602.06658. [34] T. Kerdoncuff, R. Emonet, and M. Sebban, Sampled Gromov Wasserstein, Mach. Learn., 110 (2021), pp. 2151–2186, https://doi.org/10.1007/s10994-021-06035-1. [35] B. Kloeckner, Approximation by finitely supported measures, ESAIM Control Optim. Calc. Var., 18 (2012), pp. 343–359, https://doi.org/10.1051/cocv/2010100. [36] H. Li, A. Kadav, I. Durdanovic, H. Samet, and H. P. Graf, Pruning filters for efficient ConvNets, in International Conference on Learning Representations (ICLR), 2017. [37] S. P. Lloyd, Least squares quantization in PCM, IEEE Trans. Inf. Theory, 28 (1982), pp. 129–137, https://doi.org/10.1109/TIT.1982.1056489. [38] A. Loukas, Graph reduction with spectral and cut guarantees, J. Mach. Learn. Res., 20 (2019), pp. 1–42, https://jmlr.org/papers/v20/18-680.html. [39] M. Mahajan, P. Nimbhorkar, and K. Varadarajan, The planar k-means problem is NP-hard, Theor. Comput. Sci., 442 (2012), pp. 13–21, https://doi.org/10.1016/j.tcs.2010.05.034. [40] A. W. Marshall, I. Olkin, and B. C. Arnold, Inequalities: Theory of Majorization and Its Applications, Springer Series in Statistics, Springer, New York, 2nd ed., 2011, https://doi.org/10.1007/ 978-0-387-68276-1. [41] F. Mémoli, Gromov–Wasserstein distances and the metric approach to object matching, Found. Comput. Math., 11 (2011), pp. 417–487, https://doi.org/10.1007/s10208-011-9093-5. [42] F. Mémoli and T. Needham, Comparison results for Gromov–Wasserstein and Gromov–Monge distances, ESAIM Control Optim. Calc. Var., 30 (2024), p. 78, https://doi.org/10.1051/cocv/2024063. [43] F. Mémoli, A. Sidiropoulos, and K. Singhal, Sketching and clustering metric measure spaces, 2018, https://arxiv.org/abs/1801.00551. arXiv:1801.00551.

GROMOV–WASSERSTEIN QUANTIZATION AND CLUSTERING

49

[44] F. Morelli and S. Eckstein, Partial fusion of neural networks: Efficient tradeoffs between ensembles and weight aggregation, in International Conference on Machine Learning (ICML), 2026, https:// arxiv.org/abs/2605.22350. arXiv:2605.22350. [45] S. Neumayer and G. Steidl, From optimal transport to discrepancy, in Handbook of Mathematical Models and Algorithms in Computer Vision and Imaging: Mathematical Imaging and Vision, K. Chen, C.-B. Schönlieb, X.-C. Tai, and L. Younes, eds., Springer, Cham, 2023, pp. 1791–1826, https://doi. org/10.1007/978-3-030-98661-2 95. [46] M. K. Ng, A note on constrained k-means algorithms, Pattern Recognit., 33 (2000), pp. 515–519, https: //doi.org/10.1016/S0031-3203(99)00057-6. [47] G. Pagès, Introduction to vector quantization and its applications for numerics, ESAIM Proc. Surv., 48 (2015), pp. 29–79, https://doi.org/10.1051/proc/201448002. [48] G. Peyré, M. Cuturi, and J. Solomon, Gromov–Wasserstein averaging of kernel and distance matrices, in Proceedings of the 33rd International Conference on Machine Learning, vol. 48 of Proceedings of Machine Learning Research, PMLR, 2016, pp. 2664–2672. [49] M. Piening and R. Beinert, A novel sliced fused Gromov–Wasserstein distance, in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 40, 2026, pp. 24828–24836, https://doi.org/10.1609/ aaai.v40i29.39669. [50] M. Piening and R. Beinert, Slicing Wasserstein over Wasserstein via functional optimal transport, in International Conference on Learning Representations (ICLR), 2026, https://openreview.net/forum? id=l3KtyVZde3, https://arxiv.org/abs/2509.22138. [51] L. Portales, E. Cazelles, and E. Pauwels, On the sequential convergence of Lloyd’s algorithms, Math. Oper. Res., 51 (2026), pp. 1120–1138, https://doi.org/10.1287/moor.2024.0550. [52] G. Rioux, Z. Goldfeld, and K. Kato, Semi-discrete Gromov–Wasserstein distances: Existence of Gromov–Monge maps and statistical theory, in NeurIPS 2023 Workshop on Optimal Transport and Machine Learning (OTML), 2023, https://openreview.net/forum?id=mCNl1EMAfH. [53] M. Ryner, J. Kronqvist, and J. Karlsson, Globally solving the Gromov–Wasserstein problem for point clouds in low dimensional Euclidean spaces, in Advances in Neural Information Processing Systems, vol. 36, Curran Associates, Inc., 2023, pp. 7930–7946, https://doi.org/10.52202/075280-0347. [54] M. Scetbon, G. Peyré, and M. Cuturi, Linear-time Gromov Wasserstein distances using low rank couplings and costs, in Proceedings of the 39th International Conference on Machine Learning, K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu, and S. Sabato, eds., vol. 162 of Proceedings of Machine Learning Research, PMLR, 2022, pp. 19347–19365, https://proceedings.mlr.press/v162/scetbon22b. html. [55] G. A. F. Seber and A. J. Lee, Linear Regression Analysis, Wiley Series in Probability and Statistics, Wiley, Hoboken, NJ, 2nd ed., 2003, https://doi.org/10.1002/9780471722199. [56] T. Séjourné, F.-X. Vialard, and G. Peyré, The unbalanced Gromov Wasserstein distance: Conic formulation and relaxation, in Advances in Neural Information Processing Systems, vol. 34, Curran Associates, Inc., 2021, pp. 8766–8779. [57] S. Z. Selim and M. A. Ismail, K-Means-Type algorithms: A generalized convergence theorem and characterization of local optimality, IEEE Trans. Pattern Anal. Mach. Intell., PAMI-6 (1984), pp. 81– 87, https://doi.org/10.1109/TPAMI.1984.4767478. [58] S. P. Singh and M. Jaggi, Model fusion via optimal transport, in Advances in Neural Information Processing Systems, vol. 33, Curran Associates, Inc., 2020, pp. 22045–22055. [59] J. Solomon, G. Peyré, V. G. Kim, and S. Sra, Entropic metric alignment for correspondence problems, ACM Trans. Graph., 35 (2016), pp. 72:1–72:13, https://doi.org/10.1145/2897824.2925903. [60] K.-T. Sturm, On the geometry of metric measure spaces, Acta Math., 196 (2006), pp. 65–131, https: //doi.org/10.1007/s11511-006-0002-8. [61] K.-T. Sturm, The space of spaces: curvature bounds and gradient flows on the space of metric measure spaces, Mem. Am. Math. Soc., 290 (2023), https://doi.org/10.1090/memo/1443. [62] R. W. Sumner and J. Popović, Deformation transfer for triangle meshes, ACM Trans. Graph., 23 (2004), pp. 399–405, https://doi.org/10.1145/1015706.1015736. Mesh data available at https: //people.csail.mit.edu/sumner/research/deftransfer/data.html. [63] M. Sun, Z. Liu, A. Bair, and J. Z. Kolter, A simple and effective pruning approach for large language models, in International Conference on Learning Representations (ICLR), 2024, https://openreview.

50

F. BEIER, S. ECKSTEIN

net/forum?id=PxoFut3dWW. [64] P.-N. Tan, M. Steinbach, and V. Kumar, Introduction to Data Mining, Pearson Addison-Wesley, Boston, 2006. [65] C. A. Taveras, S. Segarra, and C. A. Uribe, Gromov–Wasserstein graph coarsening, in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2026, pp. 646–650, https://doi.org/10.1109/ICASSP55912.2026.11463976, https://arxiv.org/abs/2511.08733. [66] T. Teuber, G. Steidl, P. Gwosdek, C. Schmaltz, and J. Weickert, Dithering by differences of convex functions, SIAM J. Imaging Sci., 4 (2011), pp. 79–108, https://doi.org/10.1137/100790197. [67] H. Van Assel, C. Vincent-Cuaz, N. Courty, R. Flamary, P. Frossard, and T. Vayer, Distributional reduction: Unifying dimensionality reduction and clustering with Gromov-Wasserstein, Transactions on Machine Learning Research, 2025 (2025), https://openreview.net/forum?id=cllm6SS354. [68] T. Vayer, L. Chapel, R. Flamary, R. Tavenard, and N. Courty, Fused Gromov–Wasserstein distance for structured objects, Algorithms, 13 (2020), p. 212, https://doi.org/10.3390/a13090212. [69] T. Vayer, R. Flamary, N. Courty, R. Tavenard, and L. Chapel, Sliced Gromov–Wasserstein, in Advances in Neural Information Processing Systems, vol. 32, Curran Associates, Inc., 2019, pp. 14726– 14736. [70] C. Villani, Optimal Transport: Old and New, vol. 338 of Grundlehren der mathematischen Wissenschaften, Springer, Berlin, Heidelberg, 2009, https://doi.org/10.1007/978-3-540-71050-9. [71] C. Vincent-Cuaz, R. Flamary, M. Corneli, T. Vayer, and N. Courty, Semi-relaxed Gromov– Wasserstein divergence and applications on graphs, in International Conference on Learning Representations (ICLR), 2022. [72] C. Vincent-Cuaz, T. Vayer, R. Flamary, M. Corneli, and N. Courty, Online graph dictionary learning, in Proceedings of the 38th International Conference on Machine Learning, vol. 139 of Proceedings of Machine Learning Research, PMLR, 2021, pp. 10564–10574. [73] J. von Lindheim, Simple approximative algorithms for free-support Wasserstein barycenters, Comput. Optim. Appl., 85 (2023), pp. 213–246, https://doi.org/10.1007/s10589-023-00458-3. [74] H. Xu, Gromov–Wasserstein factorization models for graph clustering, in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 34, 2020, pp. 6478–6485, https://doi.org/10.1609/aaai.v34i04. 6120. [75] H. Xu, D. Luo, and L. Carin, Scalable Gromov–Wasserstein learning for graph partitioning and matching, in Advances in Neural Information Processing Systems, vol. 32, Curran Associates, Inc., 2019, pp. 3046–3056. [76] H. Xu, D. Luo, H. Zha, and L. Carin, Gromov–Wasserstein learning for graph matching and node embedding, in Proceedings of the 36th International Conference on Machine Learning, vol. 97 of Proceedings of Machine Learning Research, PMLR, 2019, pp. 6932–6941. [77] P. L. Zador, Asymptotic quantization error of continuous signals and the quantization dimension, IEEE Trans. Inf. Theory, 28 (1982), pp. 139–149, https://doi.org/10.1109/TIT.1982.1056490. [78] Z. Zhang, Z. Goldfeld, Y. Mroueh, and B. K. Sriperumbudur, Gromov–Wasserstein distances: Entropic regularization, duality and sample complexity, Ann. Stat., 52 (2024), pp. 1616–1645, https: //doi.org/10.1214/24-AOS2406.

Related documents

• Discrete Gromov-Wasserstein Duality: Algorithms and Isomorphism Testing#925459
• Discrete Gromov-Wasserstein Duality: Algorithms and Isomorphism Testing#921223
• Averaged Mirror Descent and Dual Gradient Methods: Convergent Algorithms for Entropic Gromov-Wasserstein Problem#1115804
• Information Geometry for Wasserstein KL Divergence of Gaussian Measures on $\mathbb{R}^n$#984841
• Fast unsupervised ground metric learning with tree-Wasserstein distance#1156166