1 Introduction↩︎

In many learning tasks such as clustering or dimensionality reduction, it is often observed that the data, while embedded in a high-dimensional ambient space, actually concentrate near a low-dimensional manifold \(\mathcal{M}\) of much smaller dimension [1]–[9]. However, in practice, the geometric structure of this underlying manifold is not known in advance and has to be inferred from the data. Learning the manifold or the underlying density can be formidable challenges in high dimensions, but estimating meaningful distances that capture both geometry and density can be more tractable and sufficient for tasks like clustering. This observation has motivated extensive research on distance learning methods. Diffusion Maps [10] employ a diffusion process on the manifold and graph Laplacians to estimate diffusion distances, which capture both geometric and density information through random walks on a similarity graph. Hessian Eigenmaps [11] provide a method for recovering the underlying parametrization of scattered data lying on a manifold \(\mathcal{M}\) embedded in high-dimensional Euclidean space. Isomap [2], designed initially for dimensionality reduction, estimates geodesic distances on the manifold. Fermat distances [3] have a similar spirit but allow for metrics that depend not only on the geometry of the manifold but also on the underlying density. Another important instance where distances must be estimated arises in the Wasserstein space. In this setting, we typically do not have access to the true Wasserstein distances between probability measures but only to empirical estimates obtained from samples of the involved measures.

Despite the widespread use of these learned metrics in clustering and dimensionality reduction algorithms, the consistency of estimators constructed with them as an input has received limited attention.

The \(k\)-means clustering method provides a criterion for partitioning a cloud of points \(\{X_1, \dots, X_n \}\) into \(k\) groups: first choose cluster centers that minimize the sum of the squared distance of each point to its closest center and then assign points to each nearest cluster center. It is very fruitful to regard this procedure as a plug-in estimator of its corresponding population parameters (centers and clusters). From this perspective, Pollard, in his seminal paper [12], established the consistency of the empirical k-means for distributions on \(\mathbb{R}^d\). Later, [13] extended the consistency for Borel measures defined on separable metric spaces. The case \(k=1\) was considered previously by Ziezold in [14]. Recently, [15] introduced relaxed notions of k-means and proved consistency for their empirical counterpart for distributions defined on metric spaces that have the Heine-Borel property. In [16] the authors prove the consistency of the empirical k-means, adapted to the setting of non-uniqueness, for probability measures defined on separable Banach spaces. They also provide necessary and sufficient probabilistic conditions for the uniqueness of the k-means set of a probability distribution and determine the asymptotic distribution of the within cluster sum of squares. As a consequence, they derive a test for uniqueness of the set of k-means.

In all these references it is assumed that the metric is known in advance.

We address here the problem of stability of \(k\)-means under perturbations of the entire metric measure space. This general setting arises naturally in modern data-driven applications where the distance function is not fixed but rather evolves as a sequence that can be data-dependent (estimated through distance learning procedures) but also in randomly generated models (like first passage percolation) and discretization or approximation schemes in length spaces that can be deterministic or random.

We establish a general convergence result for both the centers and clusters obtained via the \(k\)-means procedure under such varying metrics.

Our framework considers a sequence of metric measure spaces \(\mathcal{X}_n=({\mathbb{X}}_n,d_n,\mu_n)\) and proves convergence of centers and clusters in a suitable sense, provided that the spaces \(\mathcal{X}_n\) converge to a limiting space in the measured Gromov-Hausdorff (mGH) sense, see Theorem 1 and Corollary 7.

Note that our results hold for every \(k \geq 1\) and, up to our knowledge, they are new even for \(k=1\). We then consider many applications of our main theorem to obtain the consistency of the empirical \(k\)-means and their clusters for different procedures that compute meaningful distances between the data points. This is the case of the Isomap geodesic distance estimator, Fermat distances, difussion distances arising from Laplacian eigenmaps, diffusion maps and Wasserstein barycenters computed with respect to learned ground metrics. We are not aware of any consistency results for any of them, even for \(k=1\). We also obtain new results concerning barycenters and \(k\)-barycenters in first passage percolation and for discrete approximations of length spaces.

Although our results allow for non-uniqueness of the \(k\)-means, the uniqueness problem is important for at least two reasons. On the one hand, when the set of \(k\)-means is unique, our convergence result is stronger (see equation ?? ). On the other hand, as discussed in [16], uniqueness, which is related to stability, can be used as a proxy for a good choice of \(k\). In view of this, in Subsection 3.2.3 we discuss a criterion for the uniqueness of the 1-mean under the Fermat distances. In the rest of the manuscript, we use k-means, Fréchet k-means, k-barycenters and k-medoids interchangeably.

The remainder of the paper is organized as follows. Section 2 introduces the necessary background on measured Gromov-Hausdorff convergence and presents our main consistency theorems for both centers and clusters. Section 3 contains detailed applications to the various metric learning scenarios mentioned above and develops the uniqueness criterion for Fermat distances.

2 Main results and proofs↩︎

The main result of this article is the consistency of the k-means estimator (centers and clusters) even if the notion of distance between data points is learned from the data itself rather than a given one. We remark that our statements, as well as their proofs are not data based but general abstract results on the stability of k-means under perturbations of the whole metric measure space. In section 3 we apply these results to specific tasks that involve learning from data.

We first introduce the required notions of convergence and main properties of the functionals involved in the definition of the \(k\)-means. Then we prove the convergence of the centers and the clusters under suitable mGH convergence.

2.1 Continuity of k-means with respect to the mGH topology.↩︎

We start by reviewing the measured Gromov-Hausdorff convergence (mGH). In order to do that, we first recall the notion of \(\varepsilon\)-isometry. A map \(h\colon ({\mathbb{X}},d_{\mathbb{X}}) \to ({\mathbb{Y}},d_{\mathbb{Y}})\) between metric spaces is called an \(\varepsilon\)-isometry if

• For all \(x,x' \in {\mathbb{X}}\) we have \(|d_{\mathbb{X}}(x,x') - d_{\mathbb{Y}}(h(x),h(x'))| < \varepsilon\) and

• For all \(y\in {\mathbb{Y}}\) we have \(d(y,h({\mathbb{X}}))< \varepsilon\).

A sequence of compact normalized metric measure spaces \(\mathcal{X}_n=({\mathbb{X}}_n , d_n , \mu_n )\) converges in measured Gromov–Hausdorff sense (mGH) to a compact normalized metric measure space \(\mathcal{X}=({\mathbb{X}}, d , \mu )\) if and only if there exist numbers \(\varepsilon_n \to 0\) and \(\varepsilon_n\)-isometries \(h_n\colon {\mathbb{X}}_n\to{\mathbb{X}}\) such that \((h_n)_{\#}(\mu_n) \to \mu\) weakly on \({\mathbb{X}}\) when \(n\to\infty\) ([17], [18]). We denote this convergence by \(\mathcal{X}_n \overset{mGH}{\to} \mathcal{X}\).

Note that, the existence of the \(\varepsilon_n\)-isometries \(h_n: {\mathbb{X}}_n\to {\mathbb{X}}\) is equivalent to the fact that \(\{({\mathbb{X}}_n,d_n)\}\) converges to \(({\mathbb{X}},d)\) in the Gromov-Hausdorff distance, see [19].
Given a compact metric space \((\mathbb{Y}, d_{\mathbb{Y}})\), let \(\mathcal{P}({\mathbb{Y}})\) denote the set of Borel probability measures on \(\mathbb{Y}\) equipped with the weak convergence topology. For a positive integer \(k\geq 1\), let \(C_k({\mathbb{Y}})\) denote the set of all subsets of \({\mathbb{Y}}\) of cardinality at most \(k\). For \(p\ge 1\) we define the functional \[\label{Phi} \Phi_{\mathbb{Y}}\colon C_k({\mathbb{Y}})\times \mathcal{P}({\mathbb{Y}}) \to \mathbb{R},\tag{1}\] by \[\Phi_{\mathbb{Y}}(S,\mu) = \int_{\mathbb{Y}}d_{\mathbb{Y}}^p(y,S) \, \mu(dy) = \int_{\mathbb{Y}}\min_{x\in S} d_{\mathbb{Y}}^p(y,x) \, \mu(dy).\]

For \(A\) and \(B\) in \(C_k(\mathbb{Y})\) the Hausdorff distance between them is given by \[d_H(A,B)=\max\left\{\max_{a\in A}d_{\mathbb{Y}}(a,B),\, \max_{b\in B}d_{\mathbb{Y}}(b,A)\right\},\] where \(d(x,C)=\min_{c\in C} d(x,c)\), for any \(C\in C_k({\mathbb{Y}})\). Since \({\mathbb{Y}}\) is compact, the space \((C_k(\mathbb{Y}), d_H)\) is compact as well. Lemma 4 in [9] establishes that \(\Phi_{\mathbb{Y}}\) is continuous with respect to the Hausdorff topology in the first variable and the Wasserstein topology in the second one. Thus, for each fixed \(\mu \in \mathcal{P}({\mathbb{Y}})\), \(\Phi_{\mathbb{Y}}(\cdot, \mu)\colon C_k({\mathbb{Y}}) \to \mathbb{R}\) is a continuous map on the compact space \(C_k({\mathbb{Y}})\) and therefore it attains its minimum. The elements of \(C_k({\mathbb{Y}})\) minimizing \(\Phi_{\mathbb{Y}}(\cdot, \mu)\) are called the k-means of \(\mu\). Unless strictly necessary, we do not emphasize the dependence on \(p\). As expected, the most common choice is \(p=2\).

Following [16], for a normalized metric measure space \(\mathcal{Y}=({\mathbb{Y}}, d_{\mathbb{Y}}, \mu)\) we define \[\label{k95means} \mathcal{S}_\mathcal{Y}(k):=\{S\in C_k({\mathbb{Y}}): \Phi_{\mathbb{Y}}(S,\mu)\leq \Phi_{\mathbb{Y}}(A,\mu)\; for all A\in C_k({\mathbb{Y}})\}.\tag{2}\] Observe that opposite to [16] we use the whole measure metric space \(\mathcal{Y}\) as a subscript instead of just the measure to emphasize that this functional depends also on the metric \(d_{\mathbb{Y}}\) and the space \({\mathbb{Y}}\). We can now present the main result of this work.

It is worth mentioning that the well-definiteness of both the \(\max\) and \(\min\) in ?? can be established by usual continuity and compactness arguments.

This theorem is in the same spirit as [9], with the important difference that here the distance function is unknown and is allowed to be estimated from the data. It should be interpreted as a no–false-positive guarantee: for every \(\varepsilon>0\), there exists \(N\in\mathbb{N}\) such that for all \(n\geq N\), and for every set of \(k\)-barycenters \(S_n \in \mathcal{S}_{\mathcal{X}_n}(k)\) of \(({\mathbb{X}}_n, d_n, \mu_n)\), after taking its image by the \(\varepsilon_n\)-isometry, there exists a corresponding set of \(k\)-barycenters \({S}\in \mathcal{S}_{\mathcal{X}}(k)\) of \(({\mathbb{X}}, d, \mu)\) that is \(\varepsilon\)-close in the Hausdorff distance, i.e. \[d_H(h_n(S_n), S ) \le \varepsilon.\] Before turning to the proof of Theorem 1, we state an immediate corollary, which recovers Jaffe’s results in the case of a fixed known distance [9] and will play a crucial role when applied to Wasserstein \(k\)-barycenters (see Section 3.4).

2.2 Stability of Clusters↩︎

The stability of Voronoi cells is a very natural question. However, to the best of our knowledge, it has so far been addressed only in the Euclidean setting; see Reem’s paper [22]. In [23] rates of convergence for the clusters in the \(L^2\) sense are proved for the spectral clustering algorithm on a Riemannian manifold. In this section, we prove the stability of Voronoi cells with respect to the Hausdorff topology in the general setting of compact metric spaces.

Recall that when the set of \(k\) centroids \(S=\{b_1,\cdots,\,b_k\}\) is given, the \(k\) associated clusters \(V_{b_1},\cdots, V_{b_k}\) are defined as the Voronoi cells of the centroids with respect to the set \(S\): \[\label{def:voronoicells} V_{b_i} := \left\{x\in{\mathbb{X}}\,\left|\, \forall b\in S,\, d(x,b_i)\leq d(x,b)\, \right\}\right..\tag{5}\] We use \(\mathcal{V} (S)\) for the family of Voronoi cells \(\{V_{b_1},\cdots, V_{b_k}\}\).

We can now state the main result of this section, establishing the continuity of Voronoi cells with respect to their centers.

This theorem is in the same spirit as Theorem 1. It should also be interpreted as a no–false-positive guarantee: for every \(\varepsilon>0\), there exists \(N\in\mathbb{N}\) such that for all \(n\geq N\), and for every cluster \(V\in \mathcal{V}(S_n)\), there exists a corresponding ‘true’ cluster \(W\in \mathcal{V}(S)\) whose \(\varepsilon-\)neighborhood contains the empirical cluster, i.e. \(V\subset W^\varepsilon.\)

In what follows, we present an extension of Theorem 6 to the context of clustering, contemplating Gromov-Hausdorff convergence and the existence of at most a finite numbers of elements on \(\mathcal{S}_{\mathcal{X}}(k)\).

We now turn to the proof of Theorem 6. We begin by stating the following results that will be used along its proof, included at the end of this section.

We now turn to the proof of Theorem 7. It is based on Theorem 1 which guarantees that each subsequence of \(h_n(S_n)\), with \(S_n \in \mathcal{S}_{\mathcal{X}_n}(k)\) has a further convergent subsequence whose limit point belongs to \(\mathcal{S}_{\mathcal{X}}(k)\). This can then be applied with Theorem 6 thanks to the assumption that \(\mathcal{S}_\mathcal{X}(k)\) is finite.

3 Applications↩︎

In this section, we apply our consistency results, Theorem 1 and Corollary 7, to prove that empirical k-means and their clusters do converge for some well-known metric learning procedures. We also prove convergence of \(k\)-means in other contexts with no empirical data involved like first passage percolation and discretizations of continuous length-spaces.

We start by considering the general case in which the space \({\mathbb{X}}_n\) is given by an i.i.d. sample.

3.1 Uniformly consistent distances from i.i.d. samples↩︎

The first important application of Theorem 1 concerns the general case of an i.i.d.sample \(\mathbb{X}_n = \{X_1, \dots, X_n\}\) drawn from a measure \(\mu\) on a metric space \((\mathbb{X}, d)\). More precisely, assume that \(X_i : \Omega \to \mathbb{X}\) is a sequence of i.i.d.random variables defined on a probability space \((\Omega, \mathcal{A}, P)\). By considering the empirical measure \(\mu_n\), i.e., the uniform distribution on \(\mathbb{X}_n\), we obtain a random metric measure space \((\mathbb{X}_n, d_n, \mu_n)\), where the distances \(d_n\) can be chosen in various ways. In this context, it is natural to consider, for all \(n\), the inclusion function \(h_n(x) = x\) as a family of candidates for the \(\varepsilon_n\)-isometries required for mGH convergence.

In this setting, since the empirical measures converges weakly to the population one, the mGH convergence of \((\mathbb{X}_n, d_n, \mu_n)\) towards \((\mathbb{X}, d, \mu)\) is reduced to verifying that the immersions \(h_n\colon \mathbb{X}_n \to \mathbb{X}\) give in fact a family of \(\varepsilon_n\)-isometries with \(\varepsilon_n \to 0\). That is the content of the following proposition.

3.2 Fermat and Isomap distances↩︎

Isomap [2] and Fermat distances [3], [24]–[26] have been introduced as powerful tools to estimate intrinsic distances on nonlinear data structures embedded in high-dimensional spaces.

Isomap (Isometric Mapping) was originally introduced as a nonlinear dimensionality reduction technique that aims to preserve the intrinsic geometric structure by approximating geodesic distances through neighborhood graphs. Unlike traditional linear methods such as PCA, Isomap captures the underlying manifold structure by constructing a graph where each point is connected to its nearest neighbors, with edges weighted by Euclidean distances (or by constructing any other graphs in which each node is connected to points that are close to it). The geodesic distance is then estimated as the shortest length path in this weighted graph. Typically, these distances are then used as input to multidimensional scaling (MDS) for embedding data into lower-dimensional spaces while preserving global geometric relationships.

Building upon the geodesic distance concept, Fermat distances were introduced as a generalization that considers paths minimizing a power-weighted length rather than the traditional shortest path used in Isomap and thus taking into account the intrinsic density of data points. While geodesic distances focus solely on the shortest path along the manifold, Fermat distances introduce a parameter \(\alpha \geq 1\) that controls the weighting of path lengths. For \(\alpha = 1\), the Fermat distance constructed on a graph of nearest-neighbor points coincides with the Isomap geodesic distance estimator, while for \(\alpha > 1\), it favors paths that are not necessary short in length but avoid large jumps between consecutive points of the path, thus favoring high-density regions. This makes Fermat distances particularly suitable for irregular data distributions.

It has been shown empirically that performing \(k\)-means and other clustering procedures using these data-driven distances significantly outperforms the standard \(k\)-means algorithm based on Euclidean distances [25]–[29] and gives results at the state-of-the-art of most successful methods. In particular, [27] demonstrates that, at least at the population level, this approach allows the detection of clusters with arbitrary geometries. In light of this, it is important to understand whether the centers and clusters computed at the empirical level are consistent estimators of their population counterparts.

In this subsection, we study the convergence of empirical barycenters (respectively, k-barycenters) computed with respect to the empirical Fermat distance, towards the population (respectively, \(k\)-) barycenters defined with respect to the population Fermat distance. In the same spirit, we study the barycenters computed with Isomap distances.

Let \(X_1, \dots, X_n\) be a sample of \(n\) i.i.d.points with density \(f\) with respect to the volume measure on the \(\ell\)-dimensional manifold \((\mathcal{M},g)\), i.e. \(\int_\mathcal{M}f(x)\,d\mathrm{vol}_g(x) =1,\) where in local coordinate \(d\mathrm{vol}_g(x) = \sqrt{\det g}\,dx\) stands for the volume measure on \((\mathcal{M},g)\).

We may assume that \((\mathcal{M},g)\) is embedded in an Euclidean space \(\mathbb{R}^D\) of large enough dimension.

3.2.1 Fermat distances↩︎

For \(\alpha > 1\), the empirical Fermat distance is defined as: \[d_{n,\alpha}(x,y) = \inf_\gamma \sum_{i=0}^r |X_{i+1} - X_i|^\alpha,\] where \(|\cdot |\) denotes Euclidean norm in \(\mathbb{R}^D\) and the infimum is taken over all paths \(\gamma = (X_0, \dots, X_{r+1})\) with \(X_0 = x\), \(X_{r+1} = y\).

Note that we are considering any possible path in the complete graph. Isomap distances are constructed similarly but with \(\alpha=1\). In that case it is necessary to impose to the paths that two consecutive points on them have to be nearest neighbors (or any other condition that guarantees that they are close to each other in Euclidean distance), see [2]. When \(\alpha >1\) (Fermat distances) it is not necessary to impose this restriction since optimal paths prefer to do that automatically [3].

For the population version, we consider a probability measure \(\mu\) in \((\mathcal{M},g)\) with density \(f\). The population Fermat distance \(d_{f,\alpha}\) associated to \((\mathcal{M},g, \mu)\) is the Riemannian distance induced by the conformal change of metric \(\tilde{g} = f^\kappa\, g\), where \(\kappa = (1-\alpha)/\ell.\) In other words, for \(x,y\in\mathcal{M}\) \[d_{f,\alpha}(x,y) = \inf_\gamma \int_I f^{(1-\alpha)/d}(\gamma_t) |\dot{\gamma}_t| dt.\] The infimum is taken over all piecewise smooth curves \(\gamma\colon {I=}[0, 1]\to\mathcal{M}\) with \(\gamma(0)=x\) and \(\gamma(1)=y\).

When the points \(X_i\) are i.i.d. drawn from \(f\), the empirical distances converge to the population one in the sense that there exists a constant \(\mathsf c > 0\), depending only on the parameter \(\alpha\) and the dimension \(\ell\) of \(\mathcal{M}\), such that, almost surely, \[\label{eq:cvFermat} n^{{(\alpha - 1)}/{\ell}}\, d_{n,\alpha} \underset{n \to \infty}{\longrightarrow} \mathsf c\, d_{f,\alpha}.\tag{8}\]

This convergence has been established pointwise when \(\mathcal{M}\) is isometric to the closure of an open subset of some \(\mathbb{R}^\ell\) in [3] and strengthened to uniform convergence for closed (compact and boundaryless) smooth manifolds [30].

We consider \(\mathcal{S}_{\mathcal{X}_n}(1)\), the set of Fréchet means of \(\mathcal{X}_n=(\mathbb{X}_n, d_{n,\alpha}, \mu_n)\) with respect to the empirical Fermat distance: \[\mathcal{S}_{\mathcal{X}_n}(1) = \underset{x \in \mathbb{X}_n}{\mathrm{argmin}\,} \sum_{i=1}^n d_{n,\alpha}(x, X_i)^p.\] Here \({\mathbb{X}}_n= \{X_1, \dots, X_n\}\) and we implicitly take \(\mu_n\) to be the uniform measure on \(\mathbb{X}_n\). Similarly, \(\mathcal{S}_{\mathcal{X}}(1)\) denotes the population Fréchet mean of \(\mathcal{X}\) with respect to the population Fermat distance (i.e., \(\mathcal{X}=({\mathbb{X}}, d_{f,\alpha},\mu)\) with \(\mu\) being the measure with density \(f\) on \(\mathcal{M}\)), \[\mathcal{S}_{\mathcal{X}}(1) = \underset{x \in \mathcal{M}}{\mathrm{argmin}\,}\, \int_{\mathcal{M}} d_{f,\alpha}(x, y)^p \mu(dy).\] More generally, for any \(k \in \mathbb{N}\), we consider the empirical \(k\)-means \(\mathcal{S}_{\mathcal{X}_n}(k)\), \[\mathcal{S}_{\mathcal{X}_n}(k) = \underset{S \in C_k(\mathbb{X}_n)}{\mathrm{argmin}\,} \sum_{i=1}^n d_{n,\alpha}(S, X_i)^p,\] and the population \(k\)-means, \[\mathcal{S}_{\mathcal{X}}(k) = \underset{S \in C_k(\mathcal{M})}{\mathrm{argmin}\,}\, \int_{\mathcal{M}} d_{f,\alpha}(S, y)^p \mu(dy).\]

In view of the above discussion, Theorem 1 applies, yielding the following result.

3.2.2 Isomap distances↩︎

A similar result holds for the geodesic distance estimator involved in the Isomap algorithm [31] since, as with the Fermat distances, these estimators converge uniformly.

We will focus on the \(\varepsilon\)-graph in which two points \(X_i\) and \(X_j\) are adjacent if and only if \(|X_i-X_j| \le \varepsilon\). In this graph we assign the weight \(w_{ij}=|X_i-X_j|\) to the edge \(\{X_i, X_j\}\) when \(X_i \sim X_j\).

Given \(\varepsilon>0\), the Isomap distance is defined as the \(\varepsilon\)-graph distance between \(X_i\) and \(X_j\), i.e. \[d_{{\mathbb{X}}_n,\varepsilon}(x,y) = \inf_\gamma \sum_{i=0}^r |X_{i+1}- X_i|,\] where the infimum is taken over all paths \(\gamma=(X_0, \dots, X_{r+1})\) with \(X_0=x\), \(X_{r+1}=y\), \(X_i\in \mathbb{X}_n\) for every \(1\le i \le r\) and \(|X_{i+1}-X_i|\le \varepsilon\) for every \(0\le i \le r+1\).

There is an alternative version of the Isomap distance in which instead of considering the \(\varepsilon\)-graph, the \(k\)-nearest neighbors graph (\(k\)-NN) is considered. In this graph there is an edge between \(X_i\) and \(X_j\) if one of them is among the \(k\) nearest points of the other with respect to the Euclidean distance. The following theorem can be deduced from [31]. See also [32]

Note that in [31], some additional assumptions on \(\mathcal{M}\) are made, but they are not necessary for the convergence of the estimator \(d_{{\mathbb{X}}_n, \varepsilon_n}\) to \(d_{\mathcal{M}}\).

In view of this, Proposition 13 holds also when Fermat distances are replaced with Isomap distances in the regime given in Theorem 16.

3.2.3 A uniqueness criterion for Fréchet means under the Fermat distance↩︎

As shown in [16], the uniqueness of the set of \(k\)-means is a strong indicator of an appropriate choice of \(k\). Indeed, if several distinct sets of \(k\)-means exist, this strongly suggests that \(k\) does not correspond to the true number of clusters in the data. Furthermore, uniqueness is closely related to the algorithm’s stability (see [33]). In light of this, it is important to understand the conditions under which uniqueness holds. In this section, we take a first step toward a theoretical understanding of the uniqueness issue within our metric learning framework. We introduce a criterion, based on the non-positive curvature of the Fermat metric, that guarantees the uniqueness of the Fréchet mean (i.e., the case \(k=1\) and \(p=2\) in the notation of Section 2.1) computed with respect to this metric. In other words, when this criterion is satisfied, it indicates that the data are not partitioned into several groups, but rather form a single cluster.

A metric space \((X,d)\) with non-positive curvature is one in which triangles are thinner than their Euclidean counterparts. This simple comparison property implies that, in such spaces, every probability distribution with a finite first moment admits a unique Fréchet mean, just as in the Euclidean setting, which serves as the comparison model. One of the first works linking the uniqueness of the Fréchet mean and non-positive curvature is [34]. This was later refined in [35] and [36]. Note also that such a theory has been developed for spaces with curvature bounded above, even by some positive constant, in which case an additional assumption of small support is needed to guarantee uniqueness; see [6]. More recent works have further shown that in these spaces, standard algorithms for computing Fréchet means also perform well. For instance, the empirical barycenter of a sample of a bounded random variable is sub-Gaussian, thus extending Hoeffding’s inequality to this nonlinear context (see [37] for non-positively curved spaces, and [38] for the general case of spaces with curvature bounded above). In this section, we provide a condition ensuring that a Fermat metric space is non-positively curved, and therefore enjoys barycentric properties as favorable as those of the Euclidean space.

Recall from the previous section that \(\alpha>1\), and the Fermat metric is equal to the conformal change of metric \(\tilde{g}=f^\kappa\,g\), with \(\kappa=(1-\alpha)/l<0\), \(f:\mathcal{M}\to\mathbb{R}_+\) the probability density of the model, and \(l\) the intrinsic dimension of \(\mathcal{M}\).

Note that a priori the variable \(X\) can have more than one Fréchet mean with respect to the original metric \(g\), but uniqueness appears in the Fermat metric.

As a corollary, one can deduce the uniqueness of the Fréchet mean in the Fermat distance associated to any log-concave distribution in \(\mathbb{R}^\ell\). In particular, uniqueness holds for the Gaussian distribution.

Let us conclude this section by an illustration of Corollary 18 on the Gaussian distribution. To simplify, we do the computation in dimension one. Let \(f(x)=(2\pi)^{-1/2} e^{-x^2/2}\) be the Gaussian density. The conformal Fermat metric is then given by \(\tilde{g} = f^\kappa = (2\pi)^{-\kappa/2} e^{-\kappa x^2/2}\), \(\kappa = (1-\alpha)/\ell = 1-\alpha\), since we are in dimension \(l=1\), and the Fermat distance between \(x,y\) is given by \[d_{f,\kappa}(x,y) = \int_0^1 \sqrt{\tilde{g}\left(\dot{\gamma}(t),\dot{\gamma}(t)\right)} dt = (2\pi)^{-\kappa/4}|y-x| F(x,y)\] where \(\gamma(t) = (1-t)x + ty\) and \[F(x,y) = \int_0^1 e^{-\frac{\kappa}{4}((1-t)x+ty)^2} dt.\]

So up to constant, we have that if \(Z\) is a standard normal distribution, \[\mathbb{E}\, d_{f,\kappa}(0,Z) = \int_{\mathbb{R}} (2\pi)^{-\kappa/4}|y|F(0,y) (2\pi)^{-1/2}e^{-y^2/2} dy \sim \int_\mathbb{R}|y|\int_0^1 e^{-\frac{\kappa}{4}y^2t^2} dt\,e^{-\frac{y^2}{2}}dy\] And an elementary study gives that this integral is finite if, and only if, \(\alpha\in (-1,3)\). Therefore we conclude that the normal variable \(Z\) admits a finite first moment in its own Fermat distance if, and only if the parameter \(\alpha\in[1,3)\), whence from Corollary 18 it admits a unique Fréchet mean on this range. Note that on the range \(\alpha\geq 3\), it is not integrable, and hence it does not admit a Fréchet mean.

3.3 Diffusion distances↩︎

Among nonlinear spectral dimensionality reduction techniques, Laplacian eigenmaps and diffusion maps, introduced respectively in [10], [40], are perhaps the most well-known methods. Both rely on the idea of embedding the data into a lower-dimensional Euclidean space using spectral coordinates derived from a discrete Laplace operator that encodes pairwise similarities between data points. Given a dataset \({\mathbb{X}}_n = \{X_1, \ldots, X_n\} \subset {\mathbb{X}}\subset \mathbb{R}^{\ell}\), one first defines a similarity matrix \((\eta_{ij})_{ij}\), where each entry \(\eta_{ij} = \eta(X_i, X_j)\) is computed from a positive, continuous, and symmetric kernel \(\eta : {\mathbb{X}}\times {\mathbb{X}}\to \mathbb{R}_+\). In practice, this kernel is often taken to be Gaussian, \[\label{def:Gaussiankernel} \eta(x,y) = \exp \left(-\frac{|x - y|^2}{2\sigma^2}\right).\tag{10}\] The degree matrix \(D\) is then defined as the diagonal matrix with entries \(d_i = \sum_{j=1}^n \eta_{ij}\). From these definitions, one constructs the normalized graph Laplacian \[\label{def:graphLaplace} L = I - D^{-1/2} \eta D^{-1/2}\tag{11}\] where \(I\) denotes the \(n \times n\) identity matrix and by abuse of notation we used \(\eta\) also for the matrix with entries \(\eta_{ij}\). This defines a quadratic form on \(\mathbb{R}^{n}\), given for any \(Y=(y_1, \ldots, y_n) \in \mathbb{R}^n\) by

\[Y^T L Y = \frac{1}{2} \sum_{i,j=1}^n \eta_{ij} \left( \frac{y_i}{\sqrt{d_i}} - \frac{y_j}{\sqrt{d_j}} \right)^2.\] Next, it is straightforward to verify, using the normalization by \(\sqrt{d_i}\), that the spectrum of \(L\) is contained in \([0,1]\). One can then compute the first \(k\) eigenvectors \(u_1, \ldots, u_k \in \mathbb{R}^n\) corresponding to the \(k\) smallest eigenvalues \[0 = \lambda_1 \leq \cdots \leq \lambda_k \leq 1\] of the eigenproblem \[L u = \lambda u.\] The idea is then to use these eigenvectors \[\begin{align} u_1 &= (u_1^1, \ldots, u_1^n)\in\mathbb{R}^n\\ &\vdots\\ u_k &= (u_k^1, \ldots, u_k^n)\in\mathbb{R}^n \end{align}\] to embed the dataset \({\mathbb{X}}_n\) into \(\mathbb{R}^k\). The Laplacian eigenmaps and diffusion maps methods differ in how they define the embedding. The Laplacian eigenmaps approach uses the eigenvectors directly, via the map \[\begin{align} \varphi \colon {\mathbb{X}}_n &\to \mathbb{R}^k\\ X_i &\mapsto (u_1^i, \ldots, u_k^i). \end{align}\] In contrast, diffusion maps introduce a diffusion parameter \(t>0\), defining the embedding \[\begin{align} \varphi_t : {\mathbb{X}}_n &\to \mathbb{R}^k\\ X_i &\mapsto ((1-\lambda_1)^t\, u_1^i, \ldots, (1-\lambda_k)^t\, u_k^i). \end{align}\] Note that \(\varphi_t\) converges to \(\varphi\) as \(t \to 0\). The standard spectral clustering procedure consists of applying the \(k\)-means algorithm to the embedded points \(\varphi(X_i) \in \mathbb{R}^k\), which yields \(k\) clusters \(A_1, \ldots, A_k \subset \mathbb{R}^k\). One may also work with the diffusion map \(\varphi_t\) for any \(t>0\), in which case the embedded points are \(\varphi_t(X_i)\) and the \(k\)-means algorithm produces clusters \[A_{1,t}, \ldots, A_{k,t} \subset \mathbb{R}^k.\] The corresponding clusters in the original dataset are then given by \[\varphi_t^{-1}(A_{1,t}), \;\ldots,\;\varphi_t^{-1}(A_{k,t}) \subset {\mathbb{X}}_n.\] Since \(k\)-means is performed with respect to the Euclidean distance in \(\mathbb{R}^k\), it is natural, if one wishes to perform clustering directly on the dataset \({\mathbb{X}}_n\), to consider the pull-back of the Euclidean distance via \(\varphi_t\), defined by \[\begin{align} d_t : {\mathbb{X}}_n \times {\mathbb{X}}_n &\to \mathbb{R}_+\\ (X_i, X_j) &\mapsto \left\| \varphi_t(X_i) - \varphi_t(X_j) \right\|_2. \end{align}\] Note that this only defines a pseudo-distance, since \(\varphi_t\) is usually non-injective, implying that it is possible to have \(d_t(x,y)=0\) but \(x\neq y\). This pseudo-distance \(d_t\) is however referred to as the truncated diffusion distance, and depends on \(n\). We will omit this dependence in the notation for simplicity. It satisfies the explicit formula \[\label{def:diffusiondistance} d_t(X_i, X_j)^2 = \sum_{l=1}^k (1-\lambda_l)^{2 t}\, (u_l^i - u_l^j)^2.\tag{12}\] The true (non truncated) diffusion distance corresponds to the infinite sum over all eigenvalues, and it is really a distance function.

The reader should keep in mind that in order to be fully rigorous, we have to take the quotient \({\mathbb{X}}_n/\sim\) where data points having the same image by the diffusion map \(\varphi_t\) are identified, in order for the truncated distance to be a distance function. Since this would only complicate notations without make things clearer, in the following we do not mention it anymore.

Unlike Isomap or Fermat distances, which compute a (weighted) \(\varepsilon\)-graph geodesic distance by minimizing the total weight among all paths connecting two points, the diffusion distance \(d_t\) measures the discrepancy between the distribution at time \(t\) of two random walks running on sample points \({\mathbb{X}}_n\) but started at \(X_i\) and \(X_j\) respectively. The transitions are proportional to \(\eta_{ij}\).

By construction, the diffusion map \(\varphi_t\) provides (after taking quotient) an isometry between \(({\mathbb{X}}_n,d_t)\) and \((\varphi_t({\mathbb{X}}_n),|\cdot |) \subset (\mathbb{R}^k,|\cdot | )\) that preserves the uniform measure and one can expect performing \(k\)-means in both spaces to be equivalent. Since standard spectral clustering applies the \(k\)-means algorithm to the embedded points \(\varphi_t(X_i)\) using the Euclidean distance in the ambient space \(\mathbb{R}^k\), this approach differs conceptually slightly from performing \(k\)-means in \(\mathcal{X}_n = \left( {\mathbb{X}}_n,\, d_t,\, \mu_n \right)\).

Consistency of spectral clustering has been widely studied. We return to this point by the end of the subsection.

With this perspective, we can recover the consistency of the \(k\)-means clustering procedure on \(\mathcal{X}_n\) using Theorem 1 once we prove that the sequence of random metric measure spaces \(\mathcal{X}_n\) converges almost surely in the measured Gromov–Hausdorff topology to some limiting metric measure space \(\mathcal{X}\).

For a fixed \(t \geq 0\), suppose that \({\mathbb{X}}_n\) is an i.i.d.sample drawn from a probability measure \(\mu\) supported on a closed manifold \(\mathcal{M} \subset \mathbb{R}^{\mathsf D}\).

According to Section 3.1, almost sure convergence in the measured Gromov–Hausdorff topology will hold provided that the diffusion distance \(d_t\) converges almost surely and uniformly to some distance \(\delta_t\) on \(\mathcal{M}\) as \(n \to \infty\). In other words, we seek to show the existence of a distance \(\delta_t\) on \(\mathcal{M}\) such that \[\text{almost surely,} \quad \sup_{i,j} \left| d_t(X_i, X_j) - \delta_t(X_i, X_j) \right| \underset{n \to \infty}{\longrightarrow} 0.\] From 12 , this convergence will occur if, simultaneously, the first \(k\) eigenvalues \(\lambda_1 \leq \cdots \leq \lambda_k\) of \(L\) converge, that is, \[\forall i = 1, \ldots, k, \quad \lambda_i \underset{n \to \infty}{\longrightarrow} \sigma_i,\] for some limiting sequence \(\sigma_1 \leq \cdots \leq \sigma_k\), and if the associated eigenvectors \(u_1, \ldots, u_k\) also converge. More precisely, for each \(i = 1, \ldots, k\), there should exist a function \(f_i : \mathcal{M} \to \mathbb{R}\) such that \[\text{almost surely,} \quad \sup_{j = 1, \ldots, n} \big| u_i^j - f_i(X_j) \big| \underset{n \to \infty}{\longrightarrow} 0.\] It turns out that both types of convergence have been established in the literature under mild assumptions, and we recall below a result from [41] that is sufficient for our purposes. A broader discussion of the literature on the convergence of graph Laplacians will be provided afterwards.

For the definition of compact convergence of operators, we refer the reader to [41] and the references therein. Note that the limiting operator can be expressed as \[\label{def:limitLaplacian} U(f)(x) = f(x) - \mathbb{E} \left[ \frac{\eta(x,X_1)}{\sqrt{\mathsf d(x)\, \mathsf d(X_1)}}\, f(X_1) \right]\tag{13}\] where \[\mathsf d(x) := \mathbb{E}[\eta(x,X_1)].\]

We emphasize that the simplicity assumption on the eigenvalues of \(U\) is not overly restrictive, since simplicity is a generic property: in other words, almost all operators have simple eigenvalues (consider, for instance, the case of a random matrix). For more details, we refer the reader to [42], [43].

From Theorem 19, we obtain the uniform convergence of the diffusion distances.

As stated in Section 3.1, we are therefore in a position to apply our main Theorem 1, from which we obtain the consistency of the \(k\)-means procedure for diffusion distances.

Since the consistency of the clustering procedure based on diffusion distances follows from the spectral convergence of the graph Laplacian used to define these distances, we conclude this section with a brief discussion of the literature on the convergence of graph Laplacians toward continuous operators.

Recall that the origins of spectral methods for dimensionality reduction trace back to [44], [45], where heat kernels were employed to embed a manifold into an (infinite-dimensional) Hilbert space. In practice, an effective approximation of this embedding is obtained by truncating the expansion to the first \(k\) frequencies; that is, to the first \(k\) eigenvalues, yielding a finite-dimensional representation of the manifold in \(\mathbb{R}^k\). Observe that letting \(k \to \infty\) in the definition of the diffusion distance 12 recovers precisely the original embedding of Takahashi and Bérard–Besson–Gallot.

These ideas motivate a central principle underlying spectral clustering: ideal clusters should correspond to the connected components of the manifold \(\mathcal{M}\). The eigenspace associated with the eigenvalue \(0\) of the Laplace–Beltrami operator is precisely the space of functions that are constant on each connected component. An important case is that of the Gaussian kernel defined in 10 . In that setting, the limiting operator \(U\) in 13 takes the form \[U(f)(x) = f(x) - P_{\sigma^2}(f)(x),\] where \(P_{\sigma^2}\) denotes the heat semigroup at time \(\sigma^2\) generated by the Laplace-Beltrami operator. Consequently, under a suitable choice of the bandwidth parameter \(\sigma^2 \to 0\) and an appropriate renormalization of the graph Laplacian \(L\), one indeed expects the renormalized operator to converge to the Laplace-Beltrami operator.

Let us mention that another way to define a graph Laplacian is to smooth the data using a Gaussian kernel, rather than proceeding as in 11 . For example, a normalized graph Laplacian can be defined from the data \(X_1,\ldots,X_n\) as (see e.g. [46]) \[\label{def:weigthedgraphLaplacian} \Delta_{\varepsilon,n}\, g(p) := \frac{1}{n\,\varepsilon^{d+2}} \sum_{i=1}^n \frac{\exp \left(-\frac{|X_i-p|^2}{2\varepsilon}\right)}{\sum_{j \neq i} \exp \left(-\frac{|X_i-X_j|^2}{2\varepsilon}\right)} \left( g(X_i) - g(p) \right).\tag{14}\] Such operators have the advantage of being defined directly on \(\mathcal{M}\), which makes the analysis of their convergence toward the (weighted) Laplace–Beltrami operator more natural, in contrast with the discrete operators in 11 . The convergence of those operators towards the (weighted) Laplace-Beltrami operator on \(\mathcal{M}\), as well as the convergence of the associated eigenvalues and eigenfunctions, has been extensively studied in the literature. Hein, Audibert, and von Luxburg [4] were among the first to establish the strong pointwise consistency of a family of graph Laplacians with data-dependent weights toward a weighted Laplace operator, including the case where the data lie on a submanifold. Belkin and Niyogi [47] proved the convergence in probability of all eigenvalues of the (weighted) graph Laplacian toward those of the continuous (weighted) Laplace–Beltrami operator, as well as the convergence in probability of the associated eigenfunctions in the \(L^2\) norm. Independently and almost simultaneously, Giné and Koltchinskii [46] and Hein [48] showed that, almost surely, the (weighted) graph Laplacian operator converges uniformly for a suitable class of functions toward the (weighted) Laplace–Beltrami operator. Subsequently, Singer [1] provided a convergence rate for this uniform convergence. More recently, Wormell and Reich [49] provided a refined analysis in the case of distributions supported on a hypertorus, improving upon Singer’s rate and even upgrading the \(L^2\) convergence of eigenfunctions obtained by Giné-Koltchinskii and Hein to convergence in \(L^\infty\). In parallel, García Trillos, Gerlach, Hein, and Šlepčev [23] and Calder and García Trillos [50] established quantitative rates for the Belkin-Niyogi result concerning the convergence of eigenvalues and the \(L^2\) convergence of eigenfunctions. Let us also mention the recent paper [51] by Wahl, where the study of the graph Laplacian is reformulated through kernel PCA by viewing the heat kernel of the manifold as a reproducing-kernel feature map. Even more recently, Xu and Singer [5] identified conditions ensuring pointwise convergence of the graph Laplacian in the broader framework of general metric spaces, rather than smooth manifolds.

Since in our setting we are interested in the uniform convergence of eigenfunctions, the \(L^2\)-type convergence results commonly found in the literature are not sufficient. Nevertheless, such \(L^2\) convergence can often be upgraded to uniform convergence by combining compactness arguments from functional analysis with the perturbation theory of linear operators; see, for instance, [52]. A detailed discussion of these techniques lies beyond the scope of the present work.

3.4 Wasserstein k-barycenters↩︎

Clustering probability distributions has recently received considerable attention, see for instance [53]–[55]. Most of the literature, however, has focused on computational aspects, while asymptotic guarantees have remained scarce. To the best of our knowledge, [9] provides the first published consistency result for Wasserstein \(k\)-means with \(k>1\), and does so in the general setting of metric spaces under mild assumptions. Previous works addressed only the case \(k=1\) or the Euclidean distance in \(\mathbb{R}^{D}\): the consistency of Wasserstein barycenters was established in [56] for a finite number of measures in \(\mathbb{R}^\ell\) and with great generality in [57] (infinite number of measures, even uncountable are allowed and the authors consider the Wasserstein space over a general geodesic metric space). Convergence rates were later obtained in [58] for the \(2\)-Wasserstein distance of probability measures over Euclidean space.

In this section, we address the case of the Wasserstein space built over a metric space whose distance is unknown and must therefore be learned from the data. This leads to a learned Wasserstein distance on the space of probability measures, with respect to which we compute Wasserstein \(k\)-means and establish their consistency.

Before turning to the case of an unknown distance, let us first highlight how our result naturally generalizes Jaffe’s theorem [9] in the setting of a fixed, known metric.

Assume that \(({\mathbb{X}},d)\) is a compact metric space and let \(d_W\) be the \(L^2\)-Wasserstein distance induced by the quadratic cost \(d^2\) on \({\mathbb{X}}\). The compactness of \({\mathbb{X}}\) guarantees that \((\mathcal{P}({\mathbb{X}}), d_W)\) is also a compact metric space. Interestingly, in this case \(d_W\) convergence is equivalent to weak convergence.

Consider a sequence of probabilities \((\mu_n)_{n\geq 1}\) defined on \(\mathcal{P}({\mathbb{X}})\), converging weakly towards \(\mu\) and define \[\label{spaces} \mathcal{X}_n = \big(\mathcal{P}({\mathbb{X}}),\, d_W,\, \mu_n\big), \quad \mathcal{X} = \big(\mathcal{P} ({\mathbb{X}}),\, d_W,\, \mu\big).\tag{15}\]

Since we have the weak convergence \(\mu_n \to \mu\), it is straightforward to check that \(\mathcal{X}_n\) converges in the measured Gromov-Hausdorff topology to \(\mathcal{X}\).

Therefore, Corollary 2 applies, which yields the consistency of Wasserstein \(k\)-means. In particular, this recovers the result of [56], [57] in the case \(k=1\) and Jaffe’s result [9] for all \(k \ge 1\).

To be more precise, Jaffe did not explicitly state this result in the context of Wasserstein \(k\)-barycenters, but it follows immediately from an application of his main theorem [9] to this setting.

An important and practically relevant instance of this setting is the following. We are given \(m \ge k\) probability measures \(\alpha_1,\dots,\alpha_m\) on \({\mathbb{X}}\). We wish to find a Wasserstein \(k\)-barycenter of this measures. To fit the framework presented in this section, let \(\mu\) be a discrete probability on \(\mathcal{P}({\mathbb{X}})\), assigning mass \(m^{-1}\) to \(\alpha_i\), for \(i=1, \ldots, m\). Namely, \[\label{mu95pob95wass} \mu=\frac{1}{m}\sum_{i=1}^m \delta_{\alpha_i}.\tag{16}\] We look for an element in \(\mathcal{S}_{\mathcal{X}}(k)\), for \({\mathcal{X}}\) defined in 15 , with \(\mu\) given in 16 . In practice, each distribution \(\alpha_i\) is only observed through a finite sample of \(n\) points \({\mathbb{X}}_n^i\subset{\mathbb{X}}\). A natural approach is therefore to replace each \(\alpha_i\) by its empirical measure \[\widehat \alpha_i = \frac{1}{|{\mathbb{X}}_n^i|}\sum_{x\in{\mathbb{X}}_n^i} \delta_{x},\quad for i=1, \ldots, m\] and to consider \[\mu_n = \frac{1}{m}\sum_{i=1}^m \delta_{\widehat {\alpha_i}}.\] the uniform distribution on \(\widehat \alpha_1,\dots,\widehat \alpha_m\). Since \(\widehat \alpha_i \to \alpha_i\) weakly as the sample size grows (and in \(d_W)\), we obtain the weak convergence of \(\mu_n\) to \(\mu\). Therefore, the above consistency result applies, and the empirical Wasserstein \(k\)-barycenter (or baryncenters) of \(\widehat \alpha_1,\dots,\widehat \alpha_m\) converges to the Wasserstein \(k\)-barycenter of \(\alpha_1,\ldots, \alpha_m\), in the sense of ?? .

Let us mention that we do not claim any novelty in this recovery of Jaffe’s result, since the proof of our Theorem 1 relies on [9], which establishes the continuity of the clustering functional. However, our formulation via measured Gromov-Hausdorff convergence both conceptualizes Jaffe’s result and, more importantly, provides a natural framework for a generalization to the case of unknown metrics, which is the focus of the next paragraph.
We now turn to the more challenging and realistic situation where \(d\) is not observed. In this case, we first construct a learned empirical distance \(d_n\) from the sample. Using \(d_n\), we define an empirical Wasserstein distance \(d_W^{(n)}\) on the space of probability measures over the sample \({\mathbb{X}}_n\). We then compute Wasserstein \(k\)-barycenters with respect to this learned distance. The goal is to show that this full procedure, learning a metric and then performing Wasserstein \(k\)-means, is still a consistent estimator of the population \(k\)-barycenter.

Recall that \(({\mathbb{X}},d)\) is a compact metric space, and \((\mathcal{P}({\mathbb{X}}),d_W)\) denote its associated Wasserstein space. Consider a sequence of compact metric spaces \(({\mathbb{X}}_n,d_n)\) converging to \(({\mathbb{X}},d)\) in the Gromov-Hausdorff topology. By definition, there exist numbers \(\varepsilon_n \to 0\) and \(\varepsilon_n\)-isometries \[h_n : ({\mathbb{X}}_n,d_n) \longrightarrow ({\mathbb{X}},d).\] For each \(n\), let \(d_W^{(n)}\) denote the Wasserstein distance defined between probability measured in \(({\mathbb{X}}_n, d_n)\) and consider the compact spaces \((\mathcal{P}({\mathbb{X}}_n), d_W^{(n)})\).

It is a classical result that the Gromov-Hausdorff convergence of compact metric spaces implies the Gromov-Hausdorff convergence of the associated Wasserstein spaces, see for instance [59]. We shall use the following quantitative version.

Assume in addition that we are given a probability measure \(\mu \in \mathcal{P} (\mathcal{P} ({\mathbb{X}}))\), together with a sequence \((\mu_n)_n\), with \(\mu_n\in \mathcal{P} ({\mathbb{X}}_n)\), such that \(h_n^\#\mu_n\) converges weakly to \(\mu\). Then, by Theorem 23, the metric measure spaces \[\bigl(\mathcal{P} ({\mathbb{X}}_n),\, d_W^{(n)},\, \mu_n\bigr)\] converge to \(\bigl(\mathcal{P} ({\mathbb{X}}),\, d_W,\, \mu\bigr)\) in the measured Gromov-Hausdorff topology. Consequently, Theorem 1 applies and yields the consistency of Wasserstein \(k\)-barycenters in this learned setting.
A particularly relevant situation, is when we consider \(A\), \(\mathbb{X}_n^i\), \(\widehat \alpha_i\) and \(\mu_n\) as before. Letting \({\mathbb{X}}_n = \bigcup_{i=1}^m {\mathbb{X}}_n^i\), consider any consistent estimator \(d_n\) of the ground metric \(d\) (see sections 3.2 and 3.3). We may then define \(d_W^{(n)}\) as the corresponding learned Wasserstein distance on \(\mathcal{P} ({\mathbb{X}}_n)\). Thanks to Theorem 23, this implies the following measured Gromov-Hausdorff convergence: \[\bigl(\mathcal{P} ({\mathbb{X}}_n),\, d_W^{(n)},\, \mu_n\bigr) \;\xrightarrow[n\to\infty]{mGH}\; \bigl(\mathcal{P} ({\mathbb{X}}),\, d_W,\, \mu\bigr).\] Note that in this case, the \(\varepsilon_n\)-isometries are simply the natural embeddings since \({\mathbb{X}}_n\subset{\mathbb{X}}\). Therefore, by Theorem 1, the empirical Wasserstein \(k\)-barycenter of \(\mu_n\) (computed with respect to the learned Wasserstein distance \(d_W^{(n)}\)) converges to the Wasserstein \(k\)-barycenter of \(\mu\), that is, to the barycenter of the original distributions \(\alpha_1,\dots,\alpha_m\). We conclude with the following theorem.

3.5 First passage percolation↩︎

In this subsection we prove convergence of \(k\)-barycenters for first passage Percolation (FPP) models. Up to our knowledge, these results are new even for \(k=1\).

FPP is a model proposed by Hammersley and Welsh [60] to describe the propagation of a fluid or an infection in a random media. To fix ideas, we consider the \(\ell-\)dimensional lattice \(\mathbb{L}^\ell=({\mathbb{Z}}^\ell,{\mathbb{E}}^\ell)\) in which \(x,y\in {\mathbb{Z}}^\ell\) are declared neighbors (i.e. \(\{x,y\}\in {\mathbb{E}}^\ell\)) if and only if the \(1-\)norm \(|x-y|_1=1\). We denote this by \(x\sim y\). To each edge \(e\in {\mathbb{E}}^\ell\) we assign a random passage time \(\tau_e\in\mathbb{R}_+\). The random variables \((\tau_e)_{e\in {\mathbb{E}}^\ell}\) are i.i.d. and called the passage times. For simplicity we assume here that their distribution has a density with respect to the Lebesgue measure, but this is not strictly necessary. For \(x,y \in {\mathbb{Z}}^\ell\), a path from \(x\) to \(y\) is a sequence of edges \(e_1, e_2, \dots, e_m\) in \({\mathbb{E}}^\ell\) such that consecutive edges \(e_i, e_{i+1}\) share exactly one point. The cost of such a path is \(T(\gamma):=\sum_{e\in \gamma}\tau_e\). The passage time from \(x\) to \(y\) is defined by \(T(x,y) = \inf_\gamma T(\gamma)\), the infimum is over all paths from \(x\) to \(y\). Since \({\mathbb{P}}(\tau_e=0)=0\), the pair \(({\mathbb{Z}}^\ell, T)\) is a (random) metric space. For a given Borel set \(D \subset \mathbb{R}^\ell\), we can consider FPP in \({\mathbb{X}}_n= \frac{1}{n}(n D\cap \mathbb{Z}^\ell)\) and \(d_n(x,y) = T(nx,ny)/n\), which is easily seen to be a metric space as well. We equip this metric space with \(\mu_n\), the uniform measure on \({\mathbb{X}}_n\). The translation of the celebrated limit shape theorem [61] states that \(({\mathbb{X}}_n,d_n)\) converges in Gromov-Hausdorff sense to the metric space \((D,d)\), where \(d\) is a certain (deterministic) metric that depends on the edge weights for which it is very difficult to obtain precise information. By standard weak convergence arguments, it follows that \(\mu_n\) converges weakly towards the Lebesgue measure \({\rm vol}\) on \(D\).

As a consequence, we get that the metric measure spaces \(({\mathbb{X}}_n,d_n,\mu_n)\) converge in mGH topology towards \((D,d,{\rm vol})\). It is natural to ask in this situation if the k-barycenters do converge. Since we have no information about \(d\), we cannot say much about the k-barycenter of the limit. However, from the mGH convergence and Theorem 1, we obtain that the empirical k-barycenter of the space \((\mathcal{X}_n, d_n, \mu_n)\) is a consistent estimator of it. More precisely,

Alternatively, as it is more common in FPP, we can consider for every time \(t\), the ball centered at the origin, \[B(t)=\{y \in {\mathbb{Z}}^\ell \colon T(0,y) < t\}.\] Again, for \(d_t(\frac{x}{t},\frac{y}{t})=t^{-1}T(x,y)\), we have that \(\mathcal{X}_t:=(\frac{1}{t} B(t), d_t)\) is a metric space that we also equip with the uniform measure \(\mu_t\). Note that this metric measure space admits a unique barycenter that we denote \(b_t\). The uniqueness comes from the observation that since the random passage times \(\tau_e\) admit a density with respect to the Lebesgue measure, we have that for any distinct \(x,x'\in \mathcal{X}_t\), it holds that almost surely \[\sum_{y\in B(t)} d_t^2\left(\frac{x}{t}, \frac{y}{t} \right) \neq \sum_{y\in B(t)} d_t^2\left(\frac{x'}{t},\frac{y}{t} \right).\] The convergence \(\mathcal{X}_t \underset{t\to\infty}{\longrightarrow} \mathcal{X}\) in GH sense follows from the limit shape theorem [62], where \(\mathcal{X} = (\mathcal{B}_0,d)\) for some deterministic symmetric convex compact set \(\mathcal{B}_0\in\mathbb{R}^\ell\), which induces a norm that depends on the distribution of \(\tau_e\), and \(d\) is the distance induced by this norm. In this theorem in fact Hausdorff convergence is proved for the shapes in \(\mathbb{R}^\ell\), but in our context it is equivalent to Gromov-Hausdorff convergence [63]. See also [61]. In the following theorem we extend the GH convergence to \(mGH\) convergence when both spaces are equipped with uniform measure respectively (which are different).

Here dist is the distance from a point to a set induced by \(d\) as provided by Theorem 1, but recall that \(d\) is induced by a norm in \(\mathbb{R}^\ell\) and hence this convergence to zero is equivalent to convergence to zero with the Euclidean distance from a point to a set.

3.6 Discrete approximations of length spaces↩︎

It is a well-known fact that every compact length space can be obtained as a Gromov–Hausdorff limit of finite graphs, see for instance [19]. Therefore, it is natural to look for such approximations that preserve the \(k\)-barycenter, i.e., such that the \(k\)-barycenter with respect to the (weighted) empirical measure on the approximating graph converges to the \(k\)-barycenter of the approximated length space equipped with some measure. Note that this is not trivial, since different approximating sequences \((\mathbb{X}_n, d_n)\) converging in the Gromov–Hausdorff sense to the same metric space \((\mathbb{X}, d)\) may, when equipped with the uniform measure \(\mu_n\) on \(\mathbb{X}_n\), converge to different metric measure spaces. The goal of this section is to highlight that our Theorem 1 provides a positive answer to this question, as it shows that \(k\)-barycenters are indeed preserved whenever the approximating graph is not only an approximation in the GH topology but also in the mGH topology.

We have seen in Subsection 3.1 that for any fully supported measure \(\mu\) on \(\mathbb{X}\), it is possible to obtain a metric measure space \((\mathbb{X}_n, d_n, \mu_n)\) that converges in the mGH sense towards \((\mathbb{X}, d, \mu)\) by choosing \(\mathbb{X}_n\) to be an i.i.d.sample drawn from \(\mu\), and \(d_n\) the distance induced by \(d\). However, when some information about the space \((\mathbb{X}, d)\) is available, it is preferable to choose specific (non-random) discretizations of the space, for instance, to approximate a torus with \(n\) points. A standard procedure to discretize a measure space in an optimal way is through the quantization of the measure [64]. An \(n\)th-order quantization \(q_n\) of a probability measure \(\mu\) is a probability measure supported on at most \(n\) points that solves \[\label{eq:quantization} q_n \in \underset{q \in \mathcal{P}_n}{\arg\min} \, d_W(\mu, q),\tag{17}\] where \(d_W\) denotes, the \(p-\)th order Wasserstein distance and \(\mathcal{P}_n\) is the set of probability measures on \(\mathbb{X}\) supported on at most \(n\) points. Note that this problem is equivalent to the \(k-\)means problem with \(k=n\). However, in the \(k-\)means problem, \(k\) is fixed and corresponds to the number of clusters, whereas in the quantization problem, \(k\) can vary, and in particular, the larger it is, the finer the quantization will be relative to the original measurement. It follows that if \(\mu\) is fully supported, any sequence \((q_n)_{n \ge 1}\) of minimizers in 17 satisfies \(\sup_{x \in \mathbb{X}} d(x, \mathrm{supp}\, q_n) \to 0\). This is because the minimizers do better than \(n\) i.i.d. points, for which the limit holds. As a consequence, such quantized measures can be used to construct approximations of a metric measure space in the mGH topology. More generally, let \(\varepsilon_n \to 0\), and let \(\mathbb{X}_n\) be a \(\delta_n\)-net of \((\mathbb{X}, d)\), for some \(\delta_n\) to be fixed later. We construct the \(\varepsilon_n\)-neighborhood graph \(G_n\), i.e., the graph with vertex set \(\mathbb{X}_n\), where two points \(x, y \in \mathbb{X}_n\) are connected if and only if \(d(x, y) < \varepsilon_n\). The length of each edge \(\{x, y\}\) is set to be \(d(x, y)\), and we denote by \(d_n\) the induced graph distance in \(G_n\). If \(\delta_n < \frac{1}{4} \varepsilon_n^2 / \mathrm{diam}(\mathbb{X})\), then by [19], it holds that \[\sup_{x, y \in \mathbb{X}} |d_n(x, y) - d(x, y)| \le \varepsilon_n.\] Proceeding as in Subsection 3.1, we obtain that \((\mathbb{X}_n, d_n, \mu_n) \to (\mathbb{X}, d, \mu)\) in the mGH sense if and only if \(\mu_n \to \mu\) weakly. As mentioned above, this can be achieved in several different ways, for instance, by taking \(\mu_n\) to be the uniform measure on \(\mathbb{X}_n\) when \(\mathbb{X}_n\) is chosen as an i.i.d.sample from \(\mu\). One can also take \(\mu_n\) to be the quantized measure \(q_n\), and \(\mathbb{X}_n\) to be the \(n\)-centers of this quantization \(q_n\) of \(\mu\). Note that in this latter case, where \(\mathbb{X}_n\) is chosen as the \(n\)-centers of a quantization of \(\mu\), one can instead take \(\mu_n\) to be the uniform measure on these \(n\)-centers rather than the quantized measure \(q_n\). In this situation, the sequence \((\mu_n)_{n \ge 1}\) does admit a weak limit, however, this limit is not equal to \(\mu\). We refer the reader to [64] for the Euclidean case, and to [65] for absolutely continuous measures with respect to the volume measure \({\rm vol}_g\) of a Riemannian manifold \((\mathcal{M}, g)\). To summarize, in both cases, if \(\mu\) admits a density \(\rho\) with respect to the Lebesgue measure (resp.the volume measure), then the sequence \(\bar{\mu}_n\) of probability measures that are uniform on the support of \(q_n\) converges weakly to the normalized multiple of \(\rho^{-\ell/(p+\ell)}\), where \(\ell\) denotes the dimension of the ambient space and \(p\) is the Wasserstein parameter fixed for the quantization procedure 17 . Therefore, if we want to obtain the barycenter of \(\mu\), we need to consider quantizations of the measure with density proportional to \(\rho^{-(p+\ell)/\ell}\) in order to compensate. In each of these cases, Theorem 1 applies, and we can summarize it in the following proposition.

Recall that by consistent, we always mean that for every \(\varepsilon > 0\), if \(n\) is large enough, any set of \(k\)-barycenters \(S_n\) for \(({\mathbb{X}}_n, d_n, \mu_n)\) is \(\varepsilon\)-close, in the Hausdorff distance, to some set \({S}\) of \(k\)-barycenters for \((\mathbb{X}, d, \mu)\).

References↩︎