# Matching via Distance Profiles

YoonHaeng Hur  Department of Statistics, University of Chicago  Yuehaw Khoo
Department of Statistics, University of Chicago

###### Abstract

In this paper, we introduce and study matching methods based on distance
profiles. For the matching of point clouds, the proposed method is easily
implementable by solving a linear program, circumventing the computational
obstacles of quadratic matching. Also, we propose and analyze a flexible way
to execute location-to-location matching using distance profiles. Moreover, we
provide a statistical estimation error analysis in the context of location-to-
location matching using empirical process theory. Furthermore, we apply our
method to a certain model and show its noise stability by characterizing
conditions on the noise level for the matching to be successful. Lastly, we
demonstrate the performance of the proposed method and compare it with some
existing methods using synthetic and real data.

##  1 Introduction

Matching is one of the most fundamental tasks of modern data science, which
can be easily found—under possibly different names and contexts depending on
the objects to be matched—in numerous applications. For instance, the shape
correspondence problem [Van Kaick et al., 2011] in computer graphics consists
in finding a correspondence between two geometric shapes, namely, matching the
corresponding (or related) parts of different shapes, where the shapes may be
obtained in various ways, say, by using 3D scanners or imaging equipment used
for Magnetic Resonance Imaging (MRI) or Computed Tomography (CT). Graph
matching is another example with wide-ranging applications, where the goal is
to efficiently find corresponding nodes between two (or more) graphs that may
represent structured data, such as neuronal connectivity in human brains
[Sporns et al., 2005] or social networks [Borgatti et al., 2009]. Last but not
least, unsupervised machine translation [Conneau et al., 2017] is an integral
part of modern natural language processing, where the aim is to translate one
language into another by learning a mapping that matches corresponding words.

To be more specific, let us consider the problem of matching two sets of
points, also called point clouds, say, $\\{X_{1},\ldots,X_{n}\\}$ and
$\\{Y_{1},\ldots,Y_{m}\\}$ supported on possibly different metric spaces
$(\mathcal{X},d_{\mathcal{X}})$ and $(\mathcal{Y},d_{\mathcal{Y}})$,
respectively. The most standard approach to this problem is quadratic matching
which formulates and solves a quadratic program based on pairwise distances
$\\{d_{\mathcal{X}}(X_{i},X_{k})\\}_{i,k=1}^{n}$ and
$\\{d_{\mathcal{Y}}(Y_{j},Y_{\ell})\\}_{j,\ell=1}^{m}$. First, assuming $n=m$,
suppose we want to find a one-to-one correspondence, namely, a permutation
$\pi$ on $\\{1,\ldots,n\\}$ to indicate that $X_{i}$ and $Y_{\pi(i)}$ are
matched. The idea of quadratic matching is to find $\pi$ such that
$d_{\mathcal{X}}(X_{i},X_{k})$ and $d_{\mathcal{Y}}(Y_{\pi(i)},Y_{\pi(k)})$
are as close as possible for any pair $(i,k)$, leading to the following
minimization:

$\min_{\pi\in\mathcal{S}_{n}}\sum_{i,k=1}^{n}|d_{\mathcal{X}}(X_{i},X_{k})-d_{\mathcal{Y}}(Y_{\pi(i)},Y_{\pi(k)})|^{2},$
(1.1)

where $\mathcal{S}_{n}$ is the set of all permutations on $\\{1,\ldots,n\\}$.
In words, (1.1) finds $\pi$ such that the sum of
$|d_{\mathcal{X}}(X_{i},X_{k})-d_{\mathcal{Y}}(Y_{\pi(i)},Y_{\pi(k)})|^{2}$
over all pairs $(i,k)$ is minimized. By defining matrices
$A,B\in\mathbb{R}^{n\times n}$ with entries
$A_{ik}=d_{\mathcal{X}}(X_{i},X_{k})$ and
$B_{j\ell}=d_{\mathcal{Y}}(Y_{j},Y_{\ell})$, one can deduce that (1.1) is
equivalent to
$\max_{\pi\in\mathcal{S}_{n}}\sum_{i,k=1}^{n}A_{ik}B_{\pi(i)\pi(k)}$, which is
a Quadratic Assignment Problem (QAP) [Koopmans and Beckmann, 1957, Loiola et
al., 2007]. Generally, for $n\neq m$, we often want matching that permits each
$X_{i}$ to be matched with multiple $Y_{j}$’s, which can be formalized by
couplings:

$\Pi_{n,m}:=\left\\{\gamma\in\mathbb{R}_{+}^{n\times
m}:\sum_{i=1}^{n}\gamma_{ij}=\frac{1}{m}\quad\forall
j=1,\ldots,m\quad\text{and}\quad\sum_{j=1}^{m}\gamma_{ij}=\frac{1}{n}\quad\forall
i=1,\ldots,n\right\\}.$

Here, for any element $\gamma\in\Pi_{n,m}$, which we call a coupling, its
entry $\gamma_{ij}$ can interpreted as the probability of matching $X_{i}$ and
$Y_{j}$, while the constraints require that all the $X_{i}$’s (resp.
$Y_{j}$’s) are subject to the same probability $\frac{1}{n}$ (resp.
$\frac{1}{m}$) of getting matched. Following the idea of quadratic matching
used in (1.1), we seek $\gamma$ that assigns higher probabilities to the
matching of two pairs $(X_{i},Y_{j})$ and $(X_{k},Y_{\ell})$, namely, larger
$\gamma_{ij}$ and $\gamma_{k\ell}$, if $d_{\mathcal{X}}(X_{i},X_{k})$ and
$d_{\mathcal{Y}}(Y_{j},Y_{\ell})$ are close, which can be formulated as the
following minimization:

$\min_{\gamma\in\Pi_{n,m}}\sum_{i,k=1}^{n}\sum_{j,\ell=1}^{m}|d_{\mathcal{X}}(X_{i},X_{k})-d_{\mathcal{Y}}(Y_{j},Y_{\ell})|^{2}\gamma_{ij}\gamma_{k\ell}.$
(1.2)

The objective function of (1.2), which is quadratic in the variable $\gamma$,
promotes a solution $\gamma$ with larger $\gamma_{ij}\gamma_{k\ell}$ for pairs
$(i,k)$ and $(j,\ell)$ such that
$|d_{\mathcal{X}}(X_{i},X_{k})-d_{\mathcal{Y}}(Y_{j},Y_{\ell})|$ is small.

On the theoretical side, the idea of quadratic matching has deep connections
with rich theory on isomorphisms and the Gromov-Wasserstein (GW) distance
between metric measure spaces [Mémoli, 2011], which extends the minimization
(1.2) to the setting where the point clouds are replaced by probability
measures. An increasing number of methodologies are employing such theory as a
principled approach to sophisticated tasks, ranging from classical matching
[Aflalo et al., 2015, Solomon et al., 2016, Alvarez-Melis and Jaakkola, 2018]
to statistical inference [Brécheteau, 2019, Weitkamp et al., 2022] and
generative modeling [Bunne et al., 2019, Hur et al., 2021]. On the practical
side, however, quadratic matching requires solving a quadratic program—(1.1)
or (1.2)—that is NP-hard in general. To tackle this issue, several methods
have been developed, for instance, relaxation approaches for (1.1) include
convex relaxation [Aflalo et al., 2015], semidefinite relaxation [Zhao et al.,
1998, Kezurer et al., 2015], and spectral relaxation [Leordeanu and Hebert,
2005]; for (1.2), entropic regularization [Rangarajan et al., 1999, Solomon et
al., 2016, Peyré et al., 2016] has been widely used in the literature.

Motivated by the aforementioned computational impediment to quadratic
matching, this paper introduces and studies an alternative matching method
that is intuitive and easily implementable. The idea is to associate each
point with the distance profile and match points with similar distance
profiles. More precisely, we define the distance profiles at points $X_{i}$
and $Y_{j}$ as follows:

$\hat{\mu}_{X_{i}}:=\frac{1}{n}\sum_{k=1}^{n}\delta_{d_{\mathcal{X}}(X_{i},X_{k})}\quad\text{and}\quad\hat{\nu}_{Y_{j}}:=\frac{1}{m}\sum_{\ell=1}^{m}\delta_{d_{\mathcal{Y}}(Y_{j},Y_{\ell})}.$
(1.3)

Throughout the paper, $\delta_{r}$ denotes the Dirac measure at $r$. In words,
the distance profile $\hat{\mu}_{X_{i}}$ is a probability measure on
$\mathbb{R}$ supported uniformly on the set of distances
$\\{d_{\mathcal{X}}(X_{i},X_{k})\\}_{k=1}^{n}$, thereby representing the
distribution of distances from $X_{i}$ to the points $X_{1},\ldots,X_{n}$;
analogously, the distance profile $\hat{\nu}_{Y_{j}}$ is a uniform measure on
$\mathbb{R}$ supported on the set of distances from $Y_{j}$ to the points
$Y_{1},\ldots,Y_{m}$. As mentioned earlier, our matching criterion is based on
the similarity between distance profiles, namely, we match $X_{i}$ and $Y_{j}$
if their distance profiles $\hat{\mu}_{X_{i}}$ and $\hat{\nu}_{Y_{j}}$ are
similar, where the similarity is quantified by the Wasserstein-1 distance
[Villani, 2003], that is, $W_{1}(\hat{\mu}_{X_{i}},\hat{\nu}_{Y_{j}})$. More
concretely, if we want matching represented by couplings as in (1.2), we may
find $\gamma\in\Pi_{n,m}$ such that $\gamma_{ij}$ is large if
$W_{1}(\hat{\mu}_{X_{i}},\hat{\nu}_{Y_{j}})$ small, which can be formulated as
the following minimization:

$\min_{\gamma\in\Pi_{n,m}}\sum_{i=1}^{n}\sum_{j=1}^{m}W_{1}(\hat{\mu}_{X_{i}},\hat{\nu}_{Y_{j}})\gamma_{ij}.$
(1.4)

Unlike (1.2), the minimization problem (1.4) is a linear program—which amounts
to an instance of discrete optimal transport problems [Peyré and Cuturi,
2019]—and there are several polynomial-time algorithms to solve such problems.
Similarly, for the case where $n=m$, one can deduce that the counterpart of
(1.1) can be formulated as

$\min_{\pi\in\mathcal{S}_{n}}\sum_{i=1}^{n}W_{1}(\hat{\mu}_{X_{i}},\hat{\nu}_{Y_{\pi(i)}}),$
(1.5)

which is a linear assignment problem that is solvable via polynomial-time
algorithms.

Moreover, distance profiles can be used for location-to-location matching. In
many cases, the aforementioned point clouds $X_{i}$’s and $Y_{j}$’s often
approximate—via discretization or sampling—some underlying objects represented
as abstract (possibly continuous) probability measures, say, $\mu$ and $\nu$,
supported on the metric spaces $(\mathcal{X},d_{\mathcal{X}})$ and
$(\mathcal{Y},d_{\mathcal{Y}})$, respectively. Then, the main goal is to match
these underlying objects $\mu$ and $\nu$, which requires one to decide whether
to match any query points $x$ and $y$ from the supports of $\mu$ and $\nu$,
which may be different from $X_{i}$’s and $Y_{j}$’s, respectively. In other
words, we need a principled way to decide matching between any points in the
supports of $\mu$ and $\nu$. To this end, assuming the supports of $\mu$ and
$\nu$ are $\mathcal{X}$ and $\mathcal{Y}$, respectively, we can define the
distance profiles of any $x\in\mathcal{X}$ and $y\in\mathcal{Y}$ as follows:

$\hat{\mu}_{x}:=\frac{1}{n}\sum_{k=1}^{n}\delta_{d_{\mathcal{X}}(x,X_{k})}\quad\text{and}\quad\hat{\nu}_{y}:=\frac{1}{m}\sum_{\ell=1}^{m}\delta_{d_{\mathcal{Y}}(y,Y_{\ell})}.$
(1.6)

One can see that (1.6) simply generalizes (1.3) to any $x\in\mathcal{X}$ and
$y\in\mathcal{Y}$ by considering the distributions of distances from $x$ to
$X_{i}$’s and from $y$ to $Y_{j}$’s, respectively. Then, as in the
aforementioned matching between point clouds, we can decide the matching
between $x$ and $y$ depending on how large
$W_{1}(\hat{\mu}_{x},\hat{\nu}_{y})$ is. Indeed, when the point clouds are
good approximations of the underlying objects in a sense that the associated
empirical measures $\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ and
$\hat{\nu}=\frac{1}{m}\sum_{j=1}^{m}\delta_{Y_{j}}$ are close to $\mu$ and
$\nu$, respectively, the above procedure based on
$W_{1}(\hat{\mu}_{x},\hat{\nu}_{y})$ turns out to provide a reasonable
matching between the underlying objects.

#### Our contributions

In this paper, we introduce and study matching methods based on distance
profiles. For the matching of point clouds, the proposed method is easily
implementable by solving a linear program, circumventing the computational
obstacles of quadratic matching. By formally delineating the idea of the
distance profile matching method in Section 2, we discover its interesting
connection with the optimal transport problem incorporating isomorphism-
invariance. Moreover, we propose and analyze a flexible way to execute
location-to-location matching using distance profiles in Section 3, where we
also provide a statistical estimation error analysis in the context of
location-to-location matching using empirical process theory. Furthermore, we
apply our method to a certain model and show its stability under the noise
added to the model in Section 4, characterizing conditions on the noise level
for the matching to be successful. Lastly, in Section 5, we demonstrate the
performance of the proposed method and compare it with some existing methods
using synthetic and real data.

##  2 Theoretical Foundations

###  2.1 Background

In most applications, the objects to be matched are represented as point
clouds as in the previous section. Such point clouds are often obtained by
discretizing some underlying—possibly continuous—objects (e.g., points on
smooth 3D surfaces) or by sampling from certain distributions (e.g., sets of
words from languages). Following [Mémoli, 2011], we formally conceptualize the
object, whether it is a discrete point cloud or a latent thing that is
potentially continuous, via a metric measure space.

######  Definition 1.

We call $(\mathcal{X},d_{\mathcal{X}},\mu)$ a metric measure space if
$(\mathcal{X},d_{\mathcal{X}})$ is a complete and separable metric space and
$\mu$ is a Borel probability measure on $\mathcal{X}$.111In [Mémoli, 2011],
there are additional requirements that $(\mathcal{X},d_{\mathcal{X}})$ is a
compact metric space and the support of $\mu$ must coincide with
$\mathcal{X}$. Later, we will explicitly add those requirements whenever we
need them.

As claimed in [Mémoli, 2011], metric measure spaces provide a flexible and
natural way of representing objects. For instance, in computer graphics, one
can view a smooth 3D surface as a metric measure space
$(\mathcal{X},d_{\mathcal{X}},\mu)$, where $(\mathcal{X},d_{\mathcal{X}})$ is
$\mathbb{R}^{3}$ with the standard Euclidean distance and $\mu$ is the uniform
measure supported on the surface. Meanwhile, whenever we have a point cloud
$\\{X_{1},\ldots,X_{n}\\}$ supported on a complete and separable metric space
$(\mathcal{X},d_{\mathcal{X}})$, we can view it as a metric measure space by
associating it with the empirical measure
$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$; in this case, we can
think of $(\mathcal{X},d_{\mathcal{X}},\hat{\mu})$ as a metric measure space
corresponding to $\\{X_{1},\ldots,X_{n}\\}$.

As mentioned in the previous section, the idea of quadratic matching is
related to isomorphisms between metric measure spaces, which we define as
follows based on [Sturm, 2006, Mémoli, 2011].

######  Definition 2.

Let $(\mathcal{X},d_{\mathcal{X}},\mu)$ and
$(\mathcal{Y},d_{\mathcal{Y}},\nu)$ be two metric measure spaces. We call a
measurable map $T\colon\mathcal{X}\to\mathcal{Y}$ an isomorphism (or a strong
isomorphism) if

  * (i)

$T$ is bijective,

  * (ii)

$T_{\\#}\mu=\nu$, where $T_{\\#}\mu$ is a measure on $\mathcal{Y}$ defined by
$T_{\\#}\mu(B)=\mu\\{x\in\mathcal{X}:T(x)\in B\\}$ for any Borel set
$B\subset\mathcal{Y}$.

  * (iii)

$d_{\mathcal{X}}(x,x^{\prime})=d_{\mathcal{Y}}(T(x),T(x^{\prime}))$ for any
$x,x^{\prime}\in\mathcal{X}$.

We say $(\mathcal{X},d_{\mathcal{X}},\mu)$ and
$(\mathcal{Y},d_{\mathcal{Y}},\nu)$ are isomorphic (or strongly isomorphic) if
there is an isomorphism.

To understand a connection between the isomorphism and quadratic matching,
recall from (1.1) that quadratic matching consists in finding a one-to-one
correspondence $\pi$ (permutation) between two point clouds
$\\{X_{1},\ldots,X_{n}\\}$ and $\\{Y_{1},\ldots,Y_{n}\\}$ such that
$d_{\mathcal{X}}(X_{i},X_{k})$ and $d_{\mathcal{Y}}(Y_{\pi(i)},Y_{\pi(k)})$
are as close as possible for any pair $(i,k)$. Similarly, an isomorphism $T$
between metric measure spaces $(\mathcal{X},d_{\mathcal{X}},\mu)$ and
$(\mathcal{Y},d_{\mathcal{Y}},\nu)$ represents a one-to-one correspondence
between two spaces $\mathcal{X}$ and $\mathcal{Y}$ such that the distances
$d_{\mathcal{X}}(x,x^{\prime})$ and $d_{\mathcal{Y}}(T(x),T(x^{\prime}))$ are
perfectly aligned for any $x,x^{\prime}\in\mathcal{X}$, namely, an isometry,
as implied by (i) and (iii) of Definition 2. The additional condition (ii) of
Definition 2 takes into account that matching of metric measure spaces
requires matching of the associated probability measures $\mu$ and $\nu$; it
essentially states that $T$ must preserve the infinitesimal weight
$\mu(\mathrm{d}x)$ on point $x\in\mathcal{X}$, often denoted as
$\nu(\mathrm{d}T(x))=\mu(\mathrm{d}x)$. In optimal transport theory [Villani,
2003], we say $T$ is a transport map from $\mu$ to $\nu$ if $T_{\\#}\mu=\nu$
holds. In summary, an isomorphism is a bijective transport map that is also an
isometry.

As pointed out in the previous section, the Gromov-Wasserstein distance is an
extension of (1.2) to the setting where point clouds are replaced by metric
measure spaces $(\mathcal{X},d_{\mathcal{X}},\mu)$ and
$(\mathcal{Y},d_{\mathcal{Y}},\nu)$. In such a setting, couplings are defined
as Borel probability measures on the product space
$\mathcal{X}\times\mathcal{Y}$ whose marginals on $\mathcal{X},\mathcal{Y}$
are $\mu,\nu$, respectively, which we formally define as follows.

######  Definition 3.

Let $(\mathcal{X},d_{\mathcal{X}},\mu)$ and
$(\mathcal{Y},d_{\mathcal{Y}},\nu)$ be two metric measure spaces. We call a
Borel probability measure $\gamma$ on $\mathcal{X}\times\mathcal{Y}$ a
coupling if

$\gamma(A\times\mathcal{Y})=\mu(A)\quad\forall A\overset{\text{\tiny
Borel}}{\subset}\mathcal{X}\quad\text{and}\quad\gamma(\mathcal{X}\times
B)=\nu(B)\quad\forall B\overset{\text{\tiny Borel}}{\subset}\mathcal{Y}.$
(2.1)

We denote the set of all couplings by $\Pi(\mu,\nu)$.

For a coupling $\gamma$, one may intuitively think that the infinitesimal
weight $\gamma(\mathrm{d}x\times\mathrm{d}y)$ represents the probability of
matching $x\in\mathcal{X}$ and $y\in\mathcal{Y}$. Recall from (1.2) that for a
coupling $\gamma\in\Pi_{n,m}$ in the point cloud setting, its entry
$\gamma_{ij}$ represents the probability of matching $X_{i}$ and $Y_{j}$; by
viewing the point clouds as metric measure spaces associated with the
empirical measures, one can deduce that the constraints in the definition of
$\Pi_{n,m}$ amount to the marginal constraints (2.1). Generalizing the idea of
(1.2), namely, finding a coupling that minimizes
$|d_{\mathcal{X}}(X_{i},X_{k})-d_{\mathcal{Y}}(Y_{j},Y_{\ell})|$ on average,
to the metric measure setting, one obtains the following Gromov-Wasserstein
(GW) distance [Mémoli, 2011].

######  Definition 4.

For $p\geq 1$, we define the Gromov-Wasserstein distance of order $p$ between
two metric measure spaces $(\mathcal{X},d_{\mathcal{X}},\mu)$ and
$(\mathcal{Y},d_{\mathcal{Y}},\nu)$ as follows:

$\mathrm{GW}_{p}(\mu,\nu):=\left(\inf_{\gamma\in\Pi(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{Y}}\int_{\mathcal{X}\times\mathcal{Y}}|d_{\mathcal{X}}(x,x^{\prime})-d_{\mathcal{Y}}(y,y^{\prime})|^{p}\,\mathrm{d}\gamma(x,y)\mathrm{d}\gamma(x^{\prime},y^{\prime})\right)^{1/p}.$
(2.2)

Now, one can deduce that (1.2) is equivalent to
$\mathrm{GW}_{2}^{2}(\hat{\mu},\hat{\nu})$, where $\hat{\mu}$ and $\hat{\nu}$
are the empirical measures associated with the point clouds $X_{i}$’s and
$Y_{j}$’s, respectively, manifesting clear generalization of the quadratic
matching idea to the metric measure space setting. [Mémoli, 2011] rigorously
proves that $\mathrm{GW}_{p}$ satisfies the usual metric axioms under mild
conditions. Importantly, if the supports of $\mu,\nu$ are exactly
$\mathcal{X},\mathcal{Y}$, respectively, it can be shown that
$\mathrm{GW}_{p}(\mu,\nu)=0$ holds if and only if
$(\mathcal{X},d_{\mathcal{X}},\mu)$ and $(\mathcal{Y},d_{\mathcal{Y}},\nu)$
are isomorphic.

######  Remark 1 (Optimal transport and Wasserstein distance).

The aforementioned concept of matching probability measures $\mu$ and $\nu$
via transport maps or couplings is the foundation of optimal transport theory
[Villani, 2003]. In optimal transport, the following problem—often referred to
as the Kantorovich problem—is mainly studied: given two metric measure spaces
$(\mathcal{X},d_{\mathcal{X}},\mu)$ and $(\mathcal{Y},d_{\mathcal{Y}},\nu)$,

$\inf_{\gamma\in\Pi(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{Y}}c(x,y)\,\mathrm{d}\gamma(x,y),$
(2.3)

where $c\colon\mathcal{X}\times\mathcal{Y}\to\mathbb{R}_{+}$ is a measurable
function such that $c(x,y)$ represents the unit cost for moving the mass from
$x\in\mathcal{X}$ to $y\in\mathcal{Y}$. Unlike $\mathrm{GW}_{p}$, the
objective function in (2.3) is linear in the coupling $\gamma$. In most
applications, we consider (2.3) for the case where
$(\mathcal{X},d_{\mathcal{X}})=(\mathcal{Y},d_{\mathcal{Y}})$ because letting
$c=d_{\mathcal{X}}$ yields a natural notion of cost for transport via
distances; otherwise, designing a meaningful inter-space cost function $c$ can
be nontrivial as
$(\mathcal{X},d_{\mathcal{X}})\neq(\mathcal{Y},d_{\mathcal{Y}})$. When
$(\mathcal{X},d_{\mathcal{X}})=(\mathcal{Y},d_{\mathcal{Y}})$, for $p\geq 1$,
the Wasserstein-$p$ distance between $\mu$ and $\nu$ as follows:

$W_{p}(\mu,\nu):=\left(\inf_{\gamma\in\Pi(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{X}}d_{\mathcal{X}}^{p}(x,y)\,\mathrm{d}\gamma(x,y)\right)^{1/p}.$

###  2.2 Matching via distance profiles

We formally define the distance profile given a metric measure space
$(\mathcal{X},d_{\mathcal{X}},\mu)$. For each $x\in\mathcal{X}$, we define a
probability measure on $\mathbb{R}$ that represents the distribution of the
distances from $x$.

######  Definition 5.

Let $(\mathcal{X},d_{\mathcal{X}},\mu)$ be a metric measure space. For each
$x\in\mathcal{X}$, let $\mu_{x}$ be the probability measure on $\mathbb{R}$
that represents the distribution of $d_{\mathcal{X}}(x,X)$, where $X$ is a
$\mathcal{X}$-valued random variable whose law is $\mu$. Equivalently, letting
$d_{\mathcal{X},x}\colon\mathcal{X}\to\mathbb{R}$ be the map such that
$d_{\mathcal{X},x}(\cdot)=d_{\mathcal{X}}(x,\cdot)$, we may write
$\mu_{x}=(d_{\mathcal{X},x})_{\\#}\mu$. We call $\mu_{x}$ the distance profile
of $\mu$ at point $x$.

An immediate property of distance profile is that it is isomorphism-invariant,
namely, for two metric measure spaces $(\mathcal{X},d_{\mathcal{X}},\mu)$ and
$(\mathcal{Y},d_{\mathcal{Y}},\nu)$ that are isomorphic, we have
$\mu_{x}=\nu_{T(x)}$ for any isomorphism $T\colon\mathcal{X}\to\mathcal{Y}$.

As discussed earlier, quadratic matching and the Gromov-Wasserstein distance
essentially aim to make two metric measure spaces approximately isomorphic by
aligning pairwise distances as close as possible. We use a similar approach
for distance profile matching based on the aforementioned isomorphism-
invariance, namely, given two metric measure spaces
$(\mathcal{X},d_{\mathcal{X}},\mu)$ and $(\mathcal{Y},d_{\mathcal{Y}},\nu)$,
we associate $x\in\mathcal{X}$ and $y\in\mathcal{Y}$ if the distance profiles
$\mu_{x}$ and $\nu_{y}$ are similar. In this paper, we mainly use $W_{1}$, the
Wasserstein-1 distance, to measure the similarity between probability measures
defined on $\mathbb{R}$. In summary, the basic principle of distance profile
matching is to match $x\in\mathcal{X}$ and $y\in\mathcal{Y}$ provided that
$W_{1}(\mu_{x},\nu_{y})$ is small, which can be thought of as making two
metric measure spaces approximately isomorphic by aligning the distributions
of distances. We have already seen in Section 1 how to utilize this idea to
match point clouds $\\{X_{1},\ldots,X_{n}\\}$ and $\\{Y_{1},\ldots,Y_{n}\\}$
on $(\mathcal{X},d_{\mathcal{X}})$ and $(\mathcal{Y},d_{\mathcal{Y}})$,
respectively. By letting $\hat{\mu}$ and $\hat{\nu}$ be the corresponding
empirical measures based on the point clouds, one can verify that distance
profiles $\hat{\mu}_{X_{i}}$ and $\hat{\nu}_{Y_{j}}$ under Definition 5 are
exactly (1.3). Then, we compute $W_{1}(\hat{\mu}_{X_{i}},\hat{\nu}_{Y_{j}})$
for all $i,j$ and find a coupling or a permutation that makes the distance
profiles of matched pairs similar, as formulated in (1.4) and (1.5),
respectively.

#### Optimal transport with invariance

Generalizing the idea of (1.4) to metric measure spaces
$(\mathcal{X},d_{\mathcal{X}},\mu)$ and $(\mathcal{Y},d_{\mathcal{Y}},\nu)$,
one can define the following:

$\mathcal{D}_{1}(\mu,\nu):=\inf_{\gamma\in\Pi(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{Y}}W_{1}(\mu_{x},\nu_{y})\,\mathrm{d}\gamma(x,y).$
(2.4)

First, observe that (2.4) is a natural extension of (1.4) to the metric
measure space case, where we aim to find a coupling that best aligns distance
profiles of $\mu$ and $\nu$. Interestingly, we can view (2.4) as an instance
of the Kantorovich problem (2.3), where the cost function
$c\colon\mathcal{X}\times\mathcal{Y}\to\mathbb{R}_{+}$ is given by
$c(x,y)=W_{1}(\mu_{x},\nu_{y})$, namely, the Wasserstein distance between
distance profiles at $x$ and $y$. It is worth noting that such a cost function
not only defines a meaningful inter-space cost function even when
$(\mathcal{X},d_{\mathcal{X}})\neq(\mathcal{Y},d_{\mathcal{Y}})$ but also
incorporates the isomorphism-invariance, namely,
$\mathcal{D}_{1}(U_{\\#}\mu,V_{\\#}\nu)=\mathcal{D}_{1}(\mu,\nu)$ for any
isomorphisms $U\colon\mathcal{X}\to\mathcal{X}$ and
$V\colon\mathcal{Y}\to\mathcal{Y}$ from $(\mathcal{X},d_{\mathcal{X}},\mu)$
and $(\mathcal{Y},d_{\mathcal{Y}},\nu)$ to themselves, respectively; to see
this, it suffices to observe that
$W_{1}((U_{\\#}\mu)_{U(x)},(V_{\\#}\nu)_{V(y)})=W_{1}(\mu_{x},\nu_{y})$. For
the case where $\mathcal{X}=\mathcal{Y}=\mathbb{R}^{d}$ and
$d_{\mathcal{X}}=d_{\mathcal{Y}}$ are the standard Euclidean distance,
[Alvarez-Melis et al., 2019, Grave et al., 2019] introduces the following
optimal transport formulation to enforce invariance to rigid transformations:

$\inf_{\gamma\in\Pi(\mu,\nu)}\inf_{T\in\mathcal{T}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\|T(x)-y\|_{2}^{2}\,\mathrm{d}\gamma(x,y),$
(2.5)

where $\mathcal{T}$ is a class of rigid transformations from $\mathbb{R}^{d}$
to $\mathbb{R}^{d}$, which aims at incorporating invariance to the
transformations in $\mathcal{T}$. Unlike (2.5) which requires a nontrivial
joint optimization over transformations and couplings, which can be nonconvex
in general, the right-hand side of (2.4) is a simpler linear program. Figure 1
compares the matching results by the distance profile matching and the
standard optimal transport.222The presentation of Figure 1 is inspired by
Figure 2 of [Alvarez-Melis et al., 2019]. Figure 1(a) shows the matching
obtained by solving (1.4), which leads to the correct matching of point clouds
up to a rigid transformation, whereas the standard optimal transport—namely,
replacing $W_{1}(\hat{\mu}_{X_{i}},\hat{\nu}_{Y_{j}})$ in (1.4) with
$\|X_{i}-Y_{j}\|_{2}$—fails to do so as shown in Figure 1(b).

(a) Distance Profile Matching (b) Standard Optimal Transport ($W_{1}$)

Figure 1: Matching of two point clouds, where the target points (marked with
$\times$) are obtained by applying some rigid transformation to the source
points (marked with dots); here, the true correspondence is shown by the
colors of the points. The lines of (a) and (b) show the matched pairs, shown
in green and red, respectively.

######  Remark 2.

The formulation (2.4) was first introduced by [Mémoli, 2011] as a tractable
lower bound to the Gromov-Wasserstein distance. More generally, [Mémoli, 2011]
defines the following: for $p\geq 1$,

$\mathcal{D}_{p}^{p}(\mu,\nu):=\inf_{\gamma\in\Pi(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{Y}}W_{p}^{p}(\mu_{x},\nu_{y})\,\mathrm{d}\gamma(x,y).$
(2.6)

Then, one can show that $\mathcal{D}_{p}\leq\mathrm{GW}_{p}$ always holds.
Moreover, $\mathcal{D}_{p}$ satisfies suitable metric axioms under suitable
assumptions; more precisely, it is a pseudometric as it is possible that
$\mathcal{D}_{p}(\mu,\nu)=0$ holds for two metric measure spaces that are not
isomorphic. See [Mémoli and Needham, 2022] for details.

######  Remark 3.

Lastly, we mention that one may use suitable intra-space cost functions
$c_{\mathcal{X}}\colon\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ and
$c_{\mathcal{Y}}\colon\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}$ when defining
the Gromov-Wasserstein distance $\mathrm{GW}_{p}$ and the distance profiles
$\mu_{x},\nu_{y}$. More precisely, one may replace
$d_{\mathcal{X}},d_{\mathcal{Y}}$ in (2.2) with
$c_{\mathcal{X}},c_{\mathcal{Y}}$, respectively. Similarly, one may define
$\mu_{x}=(c_{\mathcal{X},x})_{\\#}\mu$, where
$c_{\mathcal{X},x}\colon\mathcal{X}\to\mathbb{R}$ is defined by
$c_{\mathcal{X},x}(\cdot)=c_{\mathcal{X}}(x,\cdot)$. When $d_{\mathcal{X}}$ is
the standard Euclidean distance, some common choices are
$c_{\mathcal{X}}=d_{\mathcal{X}}^{2}$ and
$c_{\mathcal{X}}=\exp(-d_{\mathcal{X}}^{2})$ which amount to the squared
Euclidean distance and the heat kernel, respectively. Such variations are
useful for dealing with structured data such as networks [Chowdhury and
Mémoli, 2019]. Also, for graphs, one may encode $c_{\mathcal{X}}$ and
$c_{\mathcal{Y}}$ to denote the node connectivity; then, the resulting
distance profile $\mu_{x}=(c_{\mathcal{X},x})_{\\#}\mu$ boils down to the
degree profile used in the random graph matching literature [Ding et al.,
2021].

##  3 Location-to-Location Matching

As we briefly mentioned in Section 1, we often have point clouds
$\\{X_{1},\ldots,X_{n}\\}$ and $\\{Y_{1},\ldots,Y_{m}\\}$ that approximate
some underlying objects represented as metric measure spaces
$(\mathcal{X},d_{\mathcal{X}},\mu)$ and $(\mathcal{Y},d_{\mathcal{Y}},\nu)$,
respectively. In principle, if we had full access to the (possibly continuous)
probability measures $\mu,\nu$, we could match them via distance profiles as
explained in the previous section, say, obtain a coupling between $\mu,\nu$ by
solving (2.4). In practice, however, there are several limitations. First, it
is impossible to handle all the points from the supports of $\mu,\nu$, which
are possibly continuous spaces. Hence, finding a coupling between $\mu,\nu$
may be an exorbitant task since a coupling encodes the matching of all the
points in the supports of $\mu,\nu$. A more practical task is to decide
whether to match any query points $x$ and $y$ from the supports of $\mu$ and
$\nu$, respectively. To this end, we propose an intuitive location-to-location
matching procedure, which decides matching between $x$ and $y$ depending on
how large $W_{1}(\mu_{x},\nu_{y})$ is. A catch here is that what we have in
practice are point clouds $\\{X_{1},\ldots,X_{n}\\}$ and
$\\{Y_{1},\ldots,Y_{m}\\}$ approximating $\mu$ and $\nu$. In this case, though
we cannot compute $\mu_{x},\nu_{y}$, we can easily compute
$\hat{\mu}_{x},\hat{\nu}_{y}$ as in (1.6) and use them to decide matching
between any $x\in\mathcal{X}$ and $y\in\mathcal{Y}$, namely, match them if
$W_{1}(\hat{\mu}_{x},\hat{\nu}_{y})$ is below a certain threshold. Of course,
the validity of such a procedure depends on how well the point clouds
approximate the underlying measures $\mu,\nu$. In this section, we formally
introduce the location-to-location matching procedure and analyze the
approximation errors in terms of how close the associated empirical measures
$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ and
$\hat{\nu}=\frac{1}{m}\sum_{j=1}^{m}\delta_{Y_{j}}$ are to the underlying
measures $\mu$ and $\nu$, respectively; also, in the case where the point
clouds are i.i.d. samples from $\mu$ and $\nu$, we provide concrete upper
bounds on the errors using empirical process theory.

###  3.1 Matching principle for the isomorphic case

This section establishes the main principle for the aforementioned matching
procedure assuming the underlying metric measure spaces
$(\mathcal{X},d_{\mathcal{X}},\mu)$ and $(\mathcal{Y},d_{\mathcal{Y}},\nu)$
are isomorphic; also, we assume that there is a unique isomorphism
$T\colon\mathcal{X}\to\mathcal{Y}$. In this setting, the right matching
procedure ought to match any query points $x\in\mathcal{X}$ and
$y\in\mathcal{Y}$ such that $T(x)$ and $y$ are close, namely,
$d_{\mathcal{Y}}(T(x),y)=d_{\mathcal{X}}(x,T^{-1}(y))$ is small. For such a
pair, $W_{1}(\mu_{x},\nu_{y})$ is guaranteed to be small by the following
lemma.

######  Lemma 1.

Let $(\mathcal{X},d_{\mathcal{X}},\mu)$ be a metric measure space. Then,

$W_{1}(\mu_{x},\mu_{x^{\prime}})\leq d_{\mathcal{X}}(x,x^{\prime})\quad\forall
x,x^{\prime}\in\mathcal{X}.$

###### Proof.

Define a map $S\colon\mathcal{X}\to\mathbb{R}\times\mathbb{R}$ by letting
$S(\cdot)=(d_{\mathcal{X}}(x,\cdot),d_{\mathcal{X}}(x^{\prime},\cdot))$. Then,
one can verify that $S_{\\#}\mu$ is a coupling of $\mu_{x}$ and
$\mu_{x^{\prime}}$. Therefore,

$\begin{split}W_{1}(\mu_{x},\mu_{x^{\prime}})&=\inf_{\rho\in\Pi(\mu_{x},\mu_{x^{\prime}})}\int_{\mathbb{R}\times\mathbb{R}}|z_{1}-z_{2}|\,\mathrm{d}\rho(z_{1},z_{2})\\\
&\leq\int_{\mathbb{R}\times\mathbb{R}}|z_{1}-z_{2}|\,\mathrm{d}S_{\\#}\mu(z_{1},z_{2})\\\
&=\int_{\mathcal{X}}|d_{\mathcal{X}}(x,\bar{x})-d_{\mathcal{X}}(x^{\prime},\bar{x})|\,\mathrm{d}\mu(\bar{x})\\\
&\leq d_{\mathcal{X}}(x,x^{\prime}),\end{split}$

where the last inequality is the triangle inequality of $d_{\mathcal{X}}$. ∎

By Lemma 1, for $x\in\mathcal{X}$ and $y\in\mathcal{Y}$ such that
$d_{\mathcal{Y}}(T(x),y)=d_{\mathcal{X}}(x,T^{-1}(y))\leq\varepsilon$, we have

$W_{1}(\mu_{x},\nu_{y})=W_{1}(\mu_{x},\mu_{T^{-1}(y)})\leq
d_{\mathcal{X}}(x,T^{-1}(y))\leq\varepsilon,$

where the first equality uses $\nu_{y}=\mu_{T^{-1}(y)}$ as $T$ is an
isomorphism. The above derivation suggests the following procedure: match any
query points $x\in\mathcal{X}$ and $y\in\mathcal{Y}$ if
$W_{1}(\mu_{x},\nu_{y})\leq\varepsilon$, where $\varepsilon$ is a desired
threshold, then we can match any $x\in\mathcal{X}$ and $y\in\mathcal{Y}$ such
that $d_{\mathcal{X}}(x,T^{-1}(y))\leq\varepsilon$.

In practice, however, we usually do not have full access to $\mu,\nu$ as
mentioned earlier; instead, we have point clouds $\\{X_{1},\ldots,X_{n}\\}$
and $\\{Y_{1},\ldots,Y_{m}\\}$ such that the associated empirical measures
$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ and
$\hat{\nu}=\frac{1}{m}\sum_{j=1}^{m}\delta_{Y_{j}}$ are close to $\mu$ and
$\nu$, respectively. Therefore, we need to modify the above procedure by
computing $W_{1}(\hat{\mu}_{x},\hat{\nu}_{y})$ instead of
$W_{1}(\mu_{x},\nu_{y})$. Then, we inevitably make the following approximation
error:

$\sup_{x\in\mathcal{X},y\in\mathcal{Y}}|W_{1}(\hat{\mu}_{x},\hat{\nu}_{y})-W_{1}(\mu_{x},\nu_{y})|\leq\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})+\sup_{y\in\mathcal{Y}}W_{1}(\nu_{y},\hat{\nu}_{y})=:\mathrm{err},$

where the inequality is due to the triangle inequality of $W_{1}$. Therefore,
we have

$W_{1}(\hat{\mu}_{x},\hat{\nu}_{y})\leq W_{1}(\mu_{x},\nu_{y})+\mathrm{err},$

which suggests that we should match any $x\in\mathcal{X}$ and
$y\in\mathcal{Y}$ if
$W_{1}(\hat{\mu}_{x},\hat{\nu}_{y})\leq\varepsilon+\mathrm{err}$ to ensure
matching of $(x,y)$ satisfying
$d_{\mathcal{Y}}(T(x),y)=d_{\mathcal{X}}(x,T^{-1}(y))\leq\varepsilon$.

######  Remark 4 (Non-uniqueness of matching and symmetry).

The matching procedure presented in this section aims at a simple yet
practical goal: we want to correctly match a pair $x\in\mathcal{X}$ and
$y\in\mathcal{Y}$ such that
$d_{\mathcal{Y}}(T(x),y)=d_{\mathcal{X}}(x,T^{-1}(y))$ is small, which is a
pair that ought to be matched given that $T\colon\mathcal{X}\to\mathcal{Y}$ is
the unique isomorphism. In the case where isomorphisms are not unique, the
matching is not unique in general. For instance, suppose $\mu=\nu$ are
Gaussian, say, $\mu=\nu=N(0,I_{d})$, where
$(\mathcal{X},d_{\mathcal{X}})=(\mathcal{Y},d_{\mathcal{Y}})$ are
$\mathbb{R}^{d}$ with the standard Euclidean distance. In this case, any rigid
transformation from $\mathbb{R}^{d}$ to $\mathbb{R}^{d}$ without a translation
is an isomorphism. Particularly, we have $\mu_{x}=\mu_{x^{\prime}}$ as long as
$\|x\|_{2}=\|x^{\prime}\|_{2}$, namely, we have rotational invariance. In the
presence of such symmetry, $W_{1}(\mu_{x},\nu_{y})$ is small if
$\|x\|_{2}\approx\|y\|_{2}$, which may lead to the case where
$W_{1}(\mu_{x},\nu_{y})=0$ but $\|x-y\|_{2}$ can be arbitrarily large, say,
$x=(c,0,\ldots,0)$ and $y=(0,\ldots,0,c)$ for large $c>0$. Moreover, (2.4)
admits multiple solutions as well, i.e., optimal couplings are not unique as
well. To have unique matching, it necessarily requires more assumptions on the
metric measure spaces $(\mathcal{X},d_{\mathcal{X}},\mu)$ and
$(\mathcal{Y},d_{\mathcal{Y}},\nu)$. Later in Section 4, we analyze a certain
model where the matching is unique and study the stability of the matching
procedure under perturbations.

###  3.2 Estimation error analysis

As shown above, the following error is crucial for the aforementioned matching
procedure based on distance profiles:

$\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x}).$

This section provides upper bounds on the above error in the setting where
$\hat{\mu}$ is the empirical measure based on $X_{1},\ldots,X_{n}$ that are
i.i.d. from $\mu$. At first glance, the above error can be as large as
$W_{1}(\mu,\hat{\mu})$, namely, the approximation error of $\hat{\mu}$ under
$W_{1}$, which we summarize as the following proposition.

######  Proposition 1.

Let $(\mathcal{X},d_{\mathcal{X}},\mu)$ be a metric measure space. For any
probability measure $\hat{\mu}$ defined on $(\mathcal{X},d_{\mathcal{X}})$,

$\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})\leq W_{1}(\mu,\hat{\mu}).$

However, the above inequality is loose given that the rate of convergence of
$W_{1}(\mu_{x},\hat{\mu}_{x})$ in terms of the sample size $n$ is much faster
than that of $W_{1}(\mu,\hat{\mu})$ if $\mathcal{X}$ is a high-dimensional
space [Dudley, 1969, Fournier and Guillin, 2015]. More precisely, as
$\hat{\mu}_{x}$ can be thought of as the empirical of $\mu_{x}$ that is
defined on $\mathbb{R}$, the rate of convergence of the expectation of
$W_{1}(\mu_{x},\hat{\mu}_{x})$ is known to be $O(n^{-1/2})$.333We provide this
result in Lemma 2 for completeness; see Section A.1. On the other hand, if
$\mathcal{X}\subset\mathbb{R}^{d}$, the rate of convergence of the expectation
of $W_{1}(\mu,\hat{\mu})$ is $O(n^{-1/d})$, which can be much slower than that
of $W_{1}(\mu_{x},\hat{\mu}_{x})$ for large $d$. The main result here is that
the rate of convergence of the expectation of
$\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})$ is also $O(n^{-1/2})$,
thus is as fast as that of $W_{1}(\mu_{x},\hat{\mu}_{x})$ for fixed
$x\in\mathcal{X}$ and is much faster than that of $W_{1}(\mu,\hat{\mu})$.

######  Theorem 1.

Let $(\mathcal{X},d_{\mathcal{X}},\mu)$ be a metric measure space and
$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ be the empirical measure
based on $X_{1},\ldots,X_{n}$ that are i.i.d. from $\mu$. Suppose
$\mathcal{X}$ is a closed subset of $\mathbb{R}^{d}$ and $d_{\mathcal{X}}$ is
the standard Euclidean distance. Then,

$\mathbb{E}\left[\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})\right]\leq
C\sqrt{\frac{d+1}{n}},$

where $C$ is an absolute constant that is independent of $n$, $d$, and $\mu$.
Also, for any $\delta\in(0,1)$,

$\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})\leq
C\sqrt{\frac{d+1}{n}}+\sqrt{\frac{\log(1/\delta)}{2n}}$ (3.1)

holds with probability at least $1-\delta$.

For a general metric measure space $(\mathcal{X},d_{\mathcal{X}},\mu)$ that
may not be Euclidean, we derive upper bounds using the $\varepsilon$-covering
number $\mathcal{N}(\varepsilon,\mathcal{X},d_{\mathcal{X}})$, namely, the
smallest possible number of $\varepsilon$-balls to cover $\mathcal{X}$.

######  Theorem 2.

Let $(\mathcal{X},d_{\mathcal{X}},\mu)$ be a metric measure space and
$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ be the empirical measure
based on $X_{1},\ldots,X_{n}$ that are i.i.d. from $\mu$. Suppose
$\Delta=\sup_{x,x^{\prime}\in\mathcal{X}}d_{\mathcal{X}}(x,x^{\prime})<\infty$
and the support of $\mu$ is exactly $\mathcal{X}$. Let
$\mathcal{N}(\varepsilon,\mathcal{X},d_{\mathcal{X}})$ be the
$\varepsilon$-covering number of $\mathcal{X}$ under $d_{\mathcal{X}}$, then

$\mathbb{E}\left[\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})\right]\leq\frac{C}{\sqrt{n}}\int_{0}^{\Delta/4}\sqrt{\frac{\Delta}{\varepsilon}+\log\mathcal{N}(\varepsilon,\mathcal{X},d_{\mathcal{X}})}\,\mathrm{d}\varepsilon,$
(3.2)

where $C$ is an absolute constant that is independent of $n$, $\mathcal{X}$,
and $\mu$. Also, for any $\delta\in(0,1)$,

$\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})\leq\frac{C}{\sqrt{n}}\int_{0}^{\Delta/4}\sqrt{\frac{\Delta}{\varepsilon}+\log\mathcal{N}(\varepsilon,\mathcal{X},d_{\mathcal{X}})}\,\mathrm{d}\varepsilon+\sqrt{\frac{2\Delta^{2}\log(1/\delta)}{n}}$
(3.3)

holds with probability at least $1-\delta$.

######  Remark 5.

If one assumes that $\mathcal{X}$ is a compact subset of $\mathbb{R}^{d}$ and
$d_{\mathcal{X}}$ is the standard Euclidean distance in Theorem 2, the bounds
coincide with those of Theorem 1. To see this, it suffices to notice that
$\mathcal{N}(\varepsilon,\mathcal{X},d_{\mathcal{X}})$ scales as
$(1+2/\varepsilon)^{d}$; see Section 4.2 of [Vershynin, 2018]. The key
difference, however, is that Theorem 2 requires compactness of the support of
$\mu$, while Theorem 1 does not.

The proofs of Proposition 1, Theorem 1, and Theorem 2 are provided in Section
A.1.

######  Remark 6 (Convergence under $\mathcal{D}_{1}$).

Lastly, we briefly comment on a supplementary result regarding a different yet
related estimation error $\mathcal{D}_{1}(\mu,\hat{\mu})$, where
$\mathcal{D}_{1}$ is defined in (2.4). In general, the rate of convergence of
$\mathbb{E}\mathcal{D}_{1}(\mu,\hat{\mu})$ can be as slow as that of
$\mathbb{E}W_{1}(\mu,\hat{\mu})$. In the presence of symmetry, however, it can
be shown that the convergence of $\mathbb{E}\mathcal{D}_{1}(\mu,\hat{\mu})$ is
much faster. Importantly, if $\mu$ is rotationally invariant, we can show that
the rate of convergence of $\mathbb{E}\mathcal{D}_{1}(\mu,\hat{\mu})$ is
$O(n^{-1/2})$, thus is dimension-free. See Section A.1 for details.

##  4 Noise stability of distance profile matching

As pointed out in Remark 4, matching between metric measure spaces is not
unique in general. To have unique matching, it requires suitable assumptions
on the spaces. In this section, we apply the distance profile matching method
to a specific model such that unique matching exists and study the stability
under the noise added to it.

#### Model for unique matching

Concretely, we consider discrete distributions on the Euclidean space, which
are supported on finitely many locations, say, letting
$\mathcal{X}=\mathcal{Y}=\mathbb{R}^{d}$, we consider

$\mu=\frac{1}{n}\sum_{i=1}^{n}\delta_{\theta_{i}}\quad\text{and}\quad\nu=\frac{1}{n}\sum_{i=1}^{n}\delta_{\eta_{i}},$

where
$\theta_{1},\ldots,\theta_{n},\eta_{1},\ldots,\eta_{n}\in\mathbb{R}^{d}$. As
in Section 3.1, we consider the isomorphic case, namely, we assume
$\eta_{\pi^{\ast}(i)}=T(\theta_{i})$ for $i=1,\ldots,n$, where
$T\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ is an unknown rigid transformation
and $\pi^{\ast}$ is an unknown permutation on $\\{1,\ldots,n\\}$ representing
the correspondence, say, $\theta_{i}$ and $\eta_{\pi^{\ast}(i)}$ are to be
matched. One caveat here is that, as discussed in Remark 4, if there is some
symmetry, the matching may not be unique in general; for instance, suppose
$n=3$ and $\|\theta_{i}-\theta_{j}\|_{2}$ is the same for any $i\neq j$. To
avoid this, we assume that there is separation among the distance profiles,
where we define distance profiles using the squared Euclidean distances (cf.
Remark 3) as follows:

$\Omega:=\min_{i\neq
j}W_{1}(\mu_{i},\mu_{j}),\quad\text{where}~{}\mu_{i}=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|\theta_{i}-\theta_{k}\|_{2}^{2}}\quad\forall
i=1,\ldots,n.$ (4.1)

We assume $\Omega>0$, which not only ensures that the $n$ locations
$\theta_{1},\ldots,\theta_{n}$ in the support of $\mu$ are matched with the
$n$ locations $\eta_{\pi^{\ast}(1)},\ldots,\eta_{\pi^{\ast}(n)}$ in the
support of $\nu$, but also implies that $\pi^{\ast}$ is the unique optimal
matching in a sense that

$\pi^{\ast}=\operatorname*{arg\,min}_{\pi\in\mathcal{S}_{n}}\sum_{i=1}^{n}W_{1}(\mu_{i},\nu_{\pi(i)})=:\operatorname*{arg\,min}_{\pi\in\mathcal{S}_{n}}\mathcal{Q}(\pi),$

where $\mathcal{S}_{n}$ denotes the set of all permutations on
$\\{1,\ldots,n\\}$ and
$\nu_{i}=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|\eta_{i}-\eta_{k}\|_{2}^{2}}$. To
see this, note that $\mathcal{Q}(\pi^{\ast})=0$, while $\mathcal{Q}(\pi)>0$
for $\pi\neq\pi^{\ast}$ as $\Omega>0$.

#### Noise stability

In this model, we study the noise stability of the distance profile matching
procedure, where we now have
$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ and
$\hat{\nu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{Y_{i}}$, which are perturbed
versions of $\mu$ and $\nu$, respectively, due to the noise added to the
locations as follows:

$X_{i}=\theta_{i}+\sigma\xi_{i}\quad\text{and}\quad
Y_{i}=\eta_{i}+\sigma\zeta_{i},$

where $\sigma>0$ and $\xi_{,}\ldots,\xi_{n},\zeta_{1},\ldots,\zeta_{n}$ are
i.i.d. from $N(0,I_{d})$. For instance, for $d=3$, we can think of $X_{i}$’s
and $Y_{i}$’s as the coordinates obtained by scanning some objects of interest
but from different angles, say, front and side views, where the noise is added
to the ground truth coordinates $\theta_{i}$’s and $\eta_{i}$’s during the
scanning process. Now, given the perturbed locations $X_{i}$’s and $Y_{i}$’s,
we can still apply the distance profile matching in the following way to
obtain a permutation:

$\hat{\pi}=\operatorname*{arg\,min}_{\pi\in\mathcal{S}_{n}}\sum_{i=1}^{n}W_{1}(\hat{\mu}_{i},\hat{\nu}_{\pi(i)}),$

where $\hat{\mu}_{i}$’s and $\hat{\nu}_{i}$’s are the distance profiles based
on the noisy locations, namely,

$\hat{\mu}_{i}:=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|X_{i}-X_{k}\|_{2}^{2}}\quad\text{and}\quad\hat{\nu}_{i}:=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|Y_{i}-Y_{k}\|_{2}^{2}}.$

The main result of this section is that the above procedure is stable in the
following sense: $\hat{\pi}$ recovers $\pi^{\ast}$ with high probability if
the noise level $\sigma$ is sufficiently small compared to the separation
$\Omega$. The following theorem summarizes the above result.

######  Theorem 3.

Suppose that the locations
$\theta_{1},\ldots,\theta_{n},\eta_{1},\ldots,\eta_{n}\in\mathbb{R}^{d}$
satisfy

$\eta_{\pi^{\ast}(i)}=T(\theta_{i})\quad\forall i=1,\ldots,n,$ (4.2)

for some unknown rigid transformation $T\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$
and some unknown permutation $\pi^{\ast}$ on $\\{1,\ldots,n\\}$. Now, suppose
we have the following noisy locations:

$X_{i}=\theta_{i}+\sigma\xi_{i}\quad\text{and}\quad
Y_{i}=\eta_{i}+\sigma\zeta_{i},$ (4.3)

where $\sigma>0$ and $\xi_{1},\ldots,\xi_{n},\zeta_{1},\ldots,\zeta_{n}$ are
i.i.d. from $N(0,I_{d})$. For $i\in\\{1,\ldots,n\\}$, define the distance
profiles using the squared Euclidean distances as follows:

$\hat{\mu}_{i}:=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|X_{i}-X_{k}\|_{2}^{2}}\quad\text{and}\quad\hat{\nu}_{i}:=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|Y_{i}-Y_{k}\|_{2}^{2}}.$

Now, let

$\hat{\pi}=\operatorname*{arg\,min}_{\pi\in\mathcal{S}_{n}}\sum_{i=1}^{n}W_{1}(\hat{\mu}_{i},\hat{\nu}_{\pi(i)}),$

where $\mathcal{S}_{n}$ denotes the set of all permutations on
$\\{1,\ldots,n\\}$. Then, for any $\delta\in(0,1)$,

$\mathbb{P}(\hat{\pi}=\pi^{\ast})\geq 1-\delta$

provided

$1024\sigma^{2}\log(8n^{2}/\delta)\cdot\max\left(\Theta^{2},2\sigma^{2}d\right)\leq\Omega^{2},$

where $\Theta=\max_{k,\ell}\|\theta_{k}-\theta_{\ell}\|_{2}$ and $\Omega$ is
defined as in (4.1).

We provide the proof of Theorem 3 in Section A.2.

######  Remark 7.

For the case where there is no rigid transformation, namely, $T$ is the
identity map, the noise stability guarantee in the above model is well-studied
in the literature, known as the feature matching problem. In this case, the
simplest approach is to solve the discrete optimal transport problem—or
equivalently, a linear assignment problem—using the squared Euclidean
distances, namely,

$\hat{\pi}^{\mathrm{OT}}=\operatorname*{arg\,min}_{\pi\in\mathcal{S}_{n}}\sum_{i=1}^{n}\|X_{i}-Y_{\pi(i)}\|_{2}^{2}.$
(4.4)

It is well-studied [Collier and Dalalyan, 2016, Galstyan et al., 2022] that
such a procedure is indeed guaranteed to recover $\pi^{\ast}$ with high
probability if the separation among $\theta_{i}$’s is large enough compared to
the noise level $\sigma$. However, in the presence of an unknown rigid
transformation $T$ as we discuss in this section, such a procedure is not
guaranteed to succeed; we later demonstrate this point in Section 5.3.

##  5 Experiments

###  5.1 Location-to-location matching

We apply the location-to-location matching method discussed in Section 3 to
Gaussian mixtures. First, we consider an isomorphic case as in Section 3.1:
for simplicity, we let

$\mu=\nu=\frac{1}{K}\sum_{k=1}^{K}N(\theta_{k},\sigma^{2}I_{d}),$

where $\theta_{1},\ldots,\theta_{K}\in\mathbb{R}^{d}$ are distinct. We first
draw some query points $x_{1},\ldots,x_{J}$ and $y_{1},\ldots,y_{J}$ from
$\mu$ and $\nu$, respectively. Then, we sample $X_{1},\ldots,X_{n}$ and
$Y_{1},\ldots,Y_{n}$ from $\mu$ and $\nu$, respectively, independently of the
query points. Letting $\hat{\mu},\hat{\nu}$ be the corresponding empirical
measures, we define their distance profiles at the query points as in (1.6).
Then, for each $j\in\\{1,\ldots,J\\}=:[J]$, we match $x_{j}$ with
$y_{\pi(j)}$, where
$\pi(j)=\operatorname*{arg\,min}_{\ell\in[J]}W_{1}(\hat{\mu}_{x_{j}},\hat{\nu}_{y_{\ell}})$.
In this Gaussian mixture setting, we can evaluate the matching accuracy by
checking if matched pairs from the same mixture component; in other words, for
a matched pair $(x_{j},y_{\pi(j)})$, we say it is correctly matched if they
are from the same mixture component. Figure 2(a) shows the proportion of
correctly matched pairs, that is,

$\mathrm{accuracy}=\frac{|\\{j\in[J]:x_{j}~{}\text{and}~{}y_{\pi(j)}~{}\text{are
from the same mixture component}\\}|}{J}.$ (5.1)

As one can expect, we observe higher matching accuracy for larger sample size
because the estimation errors of $\hat{\mu},\hat{\nu}$ are smaller; meanwhile,
we can see that the matching accuracy deteriorates for large enough
$\sigma$’s, approaching $\frac{1}{K}$ which is the accuracy of randomly
matching the query points.

Next, we apply the method to a situation where the last mixture components of
$\mu$ and $\nu$ differ, namely,

$\mu=\frac{1}{K}\sum_{k=1}^{K}N(\theta_{k},\sigma^{2}I_{d})\quad\text{and}\quad\nu=\frac{1}{K}\sum_{k=1}^{K}N(\bar{\theta}_{k},\sigma^{2}I_{d}),$

where $\theta_{i}=\bar{\theta}_{i}$ for $i<k$ and
$\theta_{K}\neq\bar{\theta}_{K}$. In this case, we may ask if the query points
that are not from the outlier component $\theta_{K}$ are correctly matched. In
other words, without loss of generality, assuming for some $J^{\prime}<J$ that
$x_{1},\ldots,x_{J^{\prime}}$ are from the first $K-1$ components and
$x_{J^{\prime}+1},\ldots,x_{J}$ are from the last component, we define the
accuracy by replacing $J$ in (5.1) with $J^{\prime}$. Figure 2(b) shows the
results. Though the matching works to some extent for very small $\sigma$, one
can observe the accuracy deteriorates much faster than the isomorphic case
shown in Figure 2(a) as $\sigma$ grows.

(a) Isomorphic case (b) With an outlier

Figure 2: Location-to-location matching accuracy between Gaussian mixtures
with $d=10$ and $K=10$. (a) amounts to the isomorphic case, whereas (b)
amounts to the case where the last mixture components of $\mu$ and $\nu$
differ. For each $\sigma$, the query points are fixed ($J=100$) and the
sampling of $\hat{\mu},\hat{\nu}$ are repeated for 100 times; the lines denote
the average and the shaded regions amount to standard deviations. The dashed
vertical and horizontal lines denote
$\min_{k\neq\ell}\|\theta_{k}-\theta_{\ell}\|_{2}$ and $\frac{1}{K}$,
respectively.

###  5.2 Application to the shape correspondence problem

We apply the distance profile matching method to the shape correspondence
using CAPOD data [Papadakis, 2014]. We obtain point clouds from two 3D images
representing a human body with different poses as shown in Figure 3(a) and
Figure 3(b). Here, we color each source point, where the color is the distance
from some pivot point; letting $X_{1},\ldots,X_{n}$ be the source points and
$X_{1}$ be the pivot point (some point in the left hand colored in dark
purple), we assign a value $c_{i}=\|X_{i}-X_{1}\|_{2}$ to $X_{i}$ representing
its color. Assuming the ground truth correspondence is the identity map, say,
the target points are $Y_{1}\ldots,Y_{n}$ and $(X_{i},Y_{i})$ is the true pair
for all $i=1,\ldots,n$, we color $Y_{i}$ using $c_{i}$ as well. Now, we apply
the distance profile matching by solving (1.4). For a coupling $\gamma$,
recall that $n\gamma_{1j},\ldots,n\gamma_{nj}$ can be thought of as weights of
matching $(X_{1},Y_{j}),\ldots,(X_{n},Y_{j})$, respectively. Therefore, the
weighted average $\sum_{i=1}^{n}n\gamma_{ij}c_{i}$ of the colors encodes which
source points are matched to $Y_{j}$; the correct matching ought to give
$\gamma_{jj}=\frac{1}{n}$ and $\gamma_{ij}=0$ for $i\neq j$, thus
$\sum_{i=1}^{n}n\gamma_{ij}c_{i}=c_{j}$. Figure 3(c) and Figure 3(d) color the
target points with $\\{\sum_{i=1}^{n}n\gamma_{ij}c_{i}\\}_{j=1}^{n}$, where
$\gamma$ is a solution to the distance profile matching (1.4) and the
quadratic matching (1.2), respectively. From these, we can observe that both
methods recover the ground truth color shown in Figure 3(b) for most points
with some exceptions. Notably, points representing the head are not correctly
matched for both methods. Also, there are several mismatches for the points
representing hands; especially, quadratic matching fails to correctly match
most of the points representing left hand, whereas the proposed method matches
some of them correctly (shown by dark purple points in the left hand). To
further assess the matching result, for each $j\in\\{1,\ldots,n\\}$, we
compare the ground truth correspondence $x_{i}$ and the barycenter of the
source points matched to $y_{j}$ using the weights
$\\{n\gamma_{ij}\\}_{i=1}^{n}$ by computing the distance
$\|x_{i}-\sum_{i=1}^{n}n\gamma_{ij}x_{i}\|_{2}$. Figure 4 shows the results.
We can observe that the distance profile matching produces errors close to
$0$, while yielding relatively large errors as well. On the other hand, the
right tail of the histogram corresponding to quadratic matching is lighter and
most of the errors are distributed on $[0,0.6]$ in a more uniform fashion.

(a) Source points (b) Target points (c) Distance profile (d) Quadratic
matching

Figure 3: Point clouds representing a human body shown in the frontal view.
(a) and (b) are the source and target point clouds, respectively, representing
different poses, where the colors ($c_{i}$’s) are defined by the distance from
some point in the left hand (dark purple area); the brighter the colors are,
the farther the points are, e.g., the right foot with bright yellow colors.
(c) and (d) color the target points with the weighted average, say, color
$Y_{j}$ (ground truth color is $c_{j}$) using
$\sum_{i=1}^{n}n\gamma_{ij}c_{i}$ using $\gamma$ obtained by solving (1.4) and
(1.2), respectively. Figure 4: Histograms of the matching errors represented
as distances $\|x_{i}-\sum_{i=1}^{n}n\gamma_{ij}x_{j}\|_{2}$ using $\gamma$
obtained by solving (1.4) and (1.2) for the distance profile matching and
quadratic matching, respectively.

###  5.3 Noise stability experiments

Lastly, we provide simulation results related to the noise stability discussed
in Section 4. We first obtain $n$ locations
$\theta_{1},\ldots,\theta_{n}\in\mathbb{R}^{d}$ and define the other locations
$\eta_{1},\ldots,\eta_{n}$ as in (4.2) using some rigid transformation $T$ and
a permutation $\pi^{\ast}$. Then, we generate the points according to (4.3)
and obtain $\hat{\pi}$ following the procedure explained in Theorem 3; we
estimate $\mathbb{P}(\hat{\pi}=\pi^{\ast})$ by repeating this for 100 times.
Figure 5 shows the results, where we vary $\sigma$ and compare with the
probability $\mathbb{P}(\hat{\pi}^{\mathrm{OT}}=\pi^{\ast})$ for
$\hat{\pi}^{\mathrm{OT}}$ obtained by the procedure (4.4). As shown in Figure
5(a), when the rigid transformation is based on a rotation matrix that only
affects the first two coordinates, both methods can recover the true
permutation to some extent for small enough $\sigma$; in this case, as the
rotation matrix is relatively close to the identity matrix, the standard
optimal transport method (4.4) shows good performance. For a rotation matrix
that changes all the coordinates, however, the standard optimal transport
method fails to recover the true permutation as shown in Figure 5(b), while
the distance profile matching method performs almost as in Figure 5(a).

(a) Changing two coordinates (b) Changing all coordinates

Figure 5: Noise stability results with $d=10$ and $n=100$. The plots show the
probability $\mathbb{P}(\hat{\pi}=\pi^{\ast})$, where $\hat{\pi}$ corresponds
to the one defined in Theorem 3 for the blue lines (Distance profile) and to
(4.4) for the red lines (Standard OT); here, the probability is estimated by
repeating each setting for 100 times. For (a), the rigid transformation is
based on a rotation matrix that only affects the first two coordinates, while
(b) is based on a random rotation matrix changing all $d$ coordinates.

##  6 Discussion

In this paper, we have introduced and studied matching methods based on
distance profiles. For the matching of point clouds, the proposed method is
easily implementable by solving a linear program, circumventing the
computational obstacles of quadratic matching; also, it provides a flexible
way to execute location-to-location matching. Interestingly, the proposed
method provides an optimal transport formulation incorporating isomorphism-
invariance, which is simpler than the existing approach requiring joint
optimization. As shown by the experiments, the proposed method achieves good
empirical performance. On a theoretical side, we have provided statistical
estimation error analysis in the context of location-to-location matching
using empirical process theory; plus, we have applied our method to a certain
model and showed its stability under the noise added to the model, by
characterizing conditions on the noise level for the matching to succeed.

Lastly, we briefly comment on the limitations and possible
extensions/applications of the proposed method. First, notice that solving
(1.4) in large-scale applications, where $n,m$ are very large, can be
computationally expensive. As this problem is essentially a discrete optimal
transport problem, it would be interesting to study methods to speed up the
distance profile matching method for large-scale applications. Another
intriguing research direction is to extend the stability analysis in Section 4
to more general models that are possibly continuous, for instance, mixture
models. Lastly, the pseudometric $\mathcal{D}_{1}$ defined in (2.4) can be
used for several tasks, such as generative modeling: say, one may use
$\mathcal{D}_{1}$ as a discrepancy between the real data and the generated
data. It would be interesting to compare such an extension with recent
generative modeling methods based on the Gromov-Wasserstein distance.

## References

  * [Aflalo et al., 2015] Aflalo, Y., Bronstein, A., and Kimmel, R. (2015).  On convex relaxation of graph isomorphism.  Proceedings of the National Academy of Sciences, 112(10):2942–2947.
  * [Alvarez-Melis and Jaakkola, 2018] Alvarez-Melis, D. and Jaakkola, T. S. (2018).  Gromov-Wasserstein alignment of word embedding spaces.  arXiv preprint arXiv:1809.00013.
  * [Alvarez-Melis et al., 2019] Alvarez-Melis, D., Jegelka, S., and Jaakkola, T. S. (2019).  Towards optimal transport with global invariances.  In International Conference on Artificial Intelligence and Statistics, pages 1870–1879. PMLR.
  * [Ambrosio and Gigli, 2013] Ambrosio, L. and Gigli, N. (2013).  A User’s Guide to Optimal Transport, pages 1–155.  Springer Berlin Heidelberg.
  * [Borgatti et al., 2009] Borgatti, S. P., Mehra, A., Brass, D. J., and Labianca, G. (2009).  Network analysis in the social sciences.  Science, 323(5916):892–895.
  * [Brécheteau, 2019] Brécheteau, C. (2019).  A statistical test of isomorphism between metric-measure spaces using the distance-to-a-measure signature.  Electronic Journal of Statistics, 13(1):795 – 849.
  * [Bunne et al., 2019] Bunne, C., Alvarez-Melis, D., Krause, A., and Jegelka, S. (2019).  Learning generative models across incomparable spaces.  In International conference on machine learning.
  * [Chowdhury and Mémoli, 2019] Chowdhury, S. and Mémoli, F. (2019).  The Gromov–Wasserstein distance between networks and stable network invariants.  Information and Inference: A Journal of the IMA, 8(4):757–787.
  * [Collier and Dalalyan, 2016] Collier, O. and Dalalyan, A. S. (2016).  Minimax rates in permutation estimation for feature matching.  Journal of Machine Learning Research, 17(1):162–192.
  * [Conneau et al., 2017] Conneau, A., Lample, G., Ranzato, M., Denoyer, L., and Jégou, H. (2017).  Word translation without parallel data.  arXiv preprint arXiv:1710.04087.
  * [Ding et al., 2021] Ding, J., Ma, Z., Wu, Y., and Xu, J. (2021).  Efficient random graph matching via degree profiles.  Probability Theory and Related Fields, 179:29–115.
  * [Dudley, 1969] Dudley, R. M. (1969).  The speed of mean Glivenko-Cantelli convergence.  The Annals of Mathematical Statistics, 40(1):40–50.
  * [Dudley, 1979] Dudley, R. M. (1979).  Balls in $\mathbb{R}^{k}$ do not cut all subsets of $k+2$ points.  Advances in Mathematics, 31(3):306–308.
  * [Fournier and Guillin, 2015] Fournier, N. and Guillin, A. (2015).  On the rate of convergence in Wasserstein distance of the empirical measure.  Probability Theory and Related Fields, 162(3-4):707–738.
  * [Galstyan et al., 2022] Galstyan, T., Minasyan, A., and Dalalyan, A. S. (2022).  Optimal detection of the feature matching map in presence of noise and outliers.  Electronic Journal of Statistics, 16(2):5720–5750.
  * [Grave et al., 2019] Grave, E., Joulin, A., and Berthet, Q. (2019).  Unsupervised alignment of embeddings with Wasserstein procrustes.  In International Conference on Artificial Intelligence and Statistics, pages 1880–1890. PMLR.
  * [Hur et al., 2021] Hur, Y., Guo, W., and Liang, T. (2021).  Reversible Gromov-Monge sampler for simulation-based inference.  arXiv preprint arXiv:2109.14090.
  * [Kezurer et al., 2015] Kezurer, I., Kovalsky, S. Z., Basri, R., and Lipman, Y. (2015).  Tight relaxation of quadratic matching.  Computer Graphics Forum, 34(5):115–128.
  * [Koopmans and Beckmann, 1957] Koopmans, T. C. and Beckmann, M. (1957).  Assignment problems and the location of economic activities.  Econometrica, pages 53–76.
  * [Leordeanu and Hebert, 2005] Leordeanu, M. and Hebert, M. (2005).  A spectral technique for correspondence problems using pairwise constraints.  In International Conference on Computer Vision.
  * [Loiola et al., 2007] Loiola, E. M., De Abreu, N. M. M., Boaventura-Netto, P. O., Hahn, P., and Querido, T. (2007).  A survey for the quadratic assignment problem.  European Journal of Operational Research, 176(2):657–690.
  * [Mémoli, 2011] Mémoli, F. (2011).  Gromov–Wasserstein distances and the metric approach to object matching.  Foundations of Computational Mathematics, 11:417–487.
  * [Mémoli and Needham, 2022] Mémoli, F. and Needham, T. (2022).  Distance distributions and inverse problems for metric measure spaces.  Studies in Applied Mathematics, 149(4):943–1001.
  * [Papadakis, 2014] Papadakis, P. (2014).  The canonically posed 3d objects dataset.  In Eurographics Workshop on 3D Object Retrieval, pages 33–36.
  * [Peyré and Cuturi, 2019] Peyré, G. and Cuturi, M. (2019).  Computational optimal transport: With applications to data science.  Foundations and Trends® in Machine Learning, 11(5-6):355–607.
  * [Peyré et al., 2016] Peyré, G., Cuturi, M., and Solomon, J. (2016).  Gromov-Wasserstein averaging of kernel and distance matrices.  In International Conference on Machine Learning.
  * [Rangarajan et al., 1999] Rangarajan, A., Yuille, A., and Mjolsness, E. (1999).  Convergence properties of the softassign quadratic assignment algorithm.  Neural Computation, 11(6):1455–1474.
  * [Solomon et al., 2016] Solomon, J., Peyré, G., Kim, V. G., and Sra, S. (2016).  Entropic metric alignment for correspondence problems.  ACM Transactions on Graphics (ToG), 35(4):1–13.
  * [Sporns et al., 2005] Sporns, O., Tononi, G., and Kötter, R. (2005).  The human connectome: a structural description of the human brain.  PLoS Computational Biology, 1(4):e42.
  * [Sturm, 2006] Sturm, K.-T. (2006).  On the geometry of metric measure spaces.  Acta Mathematica, 196(1):65 – 131.
  * [Van Kaick et al., 2011] Van Kaick, O., Zhang, H., Hamarneh, G., and Cohen-Or, D. (2011).  A survey on shape correspondence.  Computer Graphics Forum, 30(6):1681–1707.
  * [Vershynin, 2018] Vershynin, R. (2018).  High-Dimensional Probability: An Introduction with Applications in Data Science.  Cambridge University Press.
  * [Villani, 2003] Villani, C. (2003).  Topics in Optimal Transportation.  American Mathematical Society.
  * [Weitkamp et al., 2022] Weitkamp, C. A., Proksch, K., Tameling, C., and Munk, A. (2022).  Distribution of distances based object matching: Asymptotic inference.  Journal of the American Statistical Association, pages 1–14.
  * [Zhao et al., 1998] Zhao, Q., Karisch, S. E., Rendl, F., and Wolkowicz, H. (1998).  Semidefinite programming relaxations for the quadratic assignment problem.  Journal of Combinatorial Optimization, 2:71–109.

##  Appendix A Supplements

###  A.1 Omitted details of Section 3.2

We first prove Proposition 1.

###### Proof of Proposition 1.

Notice that for any coupling $\gamma\in\Pi(\mu,\hat{\mu})$, we have a coupling
$D_{\\#}\gamma$ is a coupling of $\mu_{x}$ and $\hat{\mu}_{x}$, where
$D\colon\mathcal{X}\times\mathcal{X}\to\mathbb{R}\times\mathbb{R}$ defined by
$D(x_{1},x_{2})=(d_{\mathcal{X}}(x,x_{1}),d_{\mathcal{X}}(x,x_{2}))$. Then,
taking $\gamma\in\Pi(\mu,\hat{\mu})$ as an optimal coupling yielding
$W_{1}(\mu,\hat{\mu})$, we have

$\begin{split}W_{1}(\mu_{x},\hat{\mu}_{x})&=\inf_{\rho\in\Pi(\mu_{x},\hat{\mu}_{x})}\int_{\mathbb{R}\times\mathbb{R}}|z_{1}-z_{2}|\,\mathrm{d}\rho(z_{1},z_{2})\\\
&\leq\int_{\mathbb{R}\times\mathbb{R}}|z_{1}-z_{2}|\,\mathrm{d}D_{\\#}\gamma(z_{1},z_{2})\\\
&=\int_{\mathcal{X}\times\mathcal{X}}|d_{\mathcal{X}}(x,x_{1})-d_{\mathcal{X}}(x,x_{2})|\,\mathrm{d}\gamma(x_{1},x_{2})\\\
&\leq\int_{\mathcal{X}\times\mathcal{X}}d_{\mathcal{X}}(x_{1},x_{2})\,\mathrm{d}\gamma(x_{1},x_{2})\\\
&=W_{1}(\mu,\hat{\mu}),\end{split}$

where the second inequality is due to the triangle inequality of
$d_{\mathcal{X}}$. ∎

Next, we provide the proofs of Theorem 1 and Theorem 2. We first start with a
simpler case: bounds on $W_{1}(\mu_{x},\hat{\mu}_{x})$ for fixed
$x\in\mathcal{X}$.

######  Lemma 2.

Let $(\mathcal{X},d_{\mathcal{X}},\mu)$ be a metric measure space and
$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ be the empirical measure
based on $X_{1},\ldots,X_{n}$ that are i.i.d. from $\mu$. For
$x\in\mathcal{X}$,

$\mathbb{E}\left[W_{1}(\mu_{x},\hat{\mu}_{x})\right]\leq\frac{C}{\sqrt{n}},$

where $C$ is an absolute constant that is independent of $x$, $n$,
$\mathcal{X}$, and $\mu$. Also, for any $\delta\in(0,1)$,

$W_{1}(\mu_{x},\hat{\mu}_{x})\leq\frac{C}{\sqrt{n}}+\sqrt{\frac{\log(1/\delta)}{2n}}$
(A.1)

holds with probability at least $1-\delta$.

###### Proof of Lemma 2.

For $x\in\mathcal{X}$ and $r\in\mathbb{R}$, let
$\phi_{x,r}(\cdot)=1\\{d_{\mathcal{X}}(x,\cdot)\leq r\\}$, which is a function
from $\mathcal{X}$ to $\\{0,1\\}$. Let $F_{x}$ be the cumulative distribution
function of $\mu_{x}$, then $F_{x}(r)=\mu\phi_{x,r}$ by definition. Similarly,
letting $\hat{F}_{x}$ be the cumulative distribution function of
$\hat{\mu}_{x}$, we also have $\hat{F}_{x}(r)=\hat{\mu}\phi_{x,r}$. Recalling
that $W_{1}$ between two probability measures on $\mathbb{R}$ is the $L_{1}$
distance between their cumulative distribution functions, Therefore,

$W_{1}(\mu_{x},\hat{\mu}_{x})=\int_{0}^{\infty}|F_{x}(r)-\hat{F}_{x}(r)|\,\mathrm{d}r=\int_{0}^{\infty}|\mu\phi_{x,r}-\hat{\mu}\phi_{x,r}|\,\mathrm{d}r\leq\sup_{\phi\in\mathcal{F}_{x}}|\mu\phi-\hat{\mu}\phi|,$
(A.2)

where $\mathcal{F}_{x}:=\\{\phi_{x,r}:r\geq 0\\}$ for any $x\in\mathcal{X}$.
By Theorem 8.23 of [Vershynin, 2018],

$\mathbb{E}\left[W_{1}(\mu_{x},\hat{\mu}_{x})\right]\leq\mathbb{E}\left[\sup_{\phi\in\mathcal{F}_{x}}|\mu\phi-\hat{\mu}\phi|\right]\leq
C\sqrt{\frac{\mathrm{vc}(\mathcal{F}_{x})}{n}}=\frac{C}{\sqrt{n}},$ (A.3)

where $C$ is an absolute constant that is independent of $x$, $n$,
$\mathcal{X}$, and $\mu$ and $\mathrm{vc}(\mathcal{F}_{x})=1$, namely, the
Vapnik-Chervonenkis (VC) dimension of a function class $\mathcal{F}_{x}$ is
$1$ for all $x\in\mathcal{X}$ by construction. Note that
$\sup_{\phi\in\mathcal{F}_{x}}|\mu\phi-\hat{\mu}\phi|$ is a random variable
that is a function of $X_{1},\ldots,X_{n}$, say,

$L(X_{1},\ldots,X_{n}):=\sup_{\phi\in\mathcal{F}_{x}}|\mu\phi-\hat{\mu}\phi|.$

One can verify that
$L\colon\mathcal{X}\times\cdots\times\mathcal{X}\to\mathbb{R}$ satisfies the
bounded difference inequality, namely, for any
$x_{1},\ldots,x_{n},x_{i}^{\prime}\in\mathcal{X}$,

$\left|L(x_{1},\ldots,x_{i},\ldots,x_{n})-L(x_{1},\ldots,x_{i}^{\prime},\ldots,x_{n})\right|\leq\frac{1}{n}.$

Therefore, by McDiarmid’s inequality (Theorem 2.9.1 of [Vershynin, 2018]),

$\sup_{\phi\in\mathcal{F}_{x}}|\mu\phi-\hat{\mu}\phi|\leq\mathbb{E}\left[\sup_{\phi\in\mathcal{F}_{x}}|\mu\phi-\hat{\mu}\phi|\right]+\sqrt{\frac{\log(1/\delta)}{2n}}$
(A.4)

holds with probability at least $1-\delta$. We obtain (A.1) by combining (A.4)
with (A.2) and (A.3). ∎

Now, we prove Theorem 1.

###### Proof of Theorem 1.

As in the proof of Lemma 2,

$\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})\leq\sup_{x\in\mathcal{X}}\sup_{\phi\in\mathcal{F}_{x}}|\mu\phi-\hat{\mu}\phi|=\sup_{\phi\in\mathcal{F}_{\mathcal{X}}}|\mu\phi-\hat{\mu}\phi|,$
(A.5)

where
$\mathcal{F}_{\mathcal{X}}:=\bigcup_{x\in\mathcal{X}}\mathcal{F}_{x}=\\{\phi_{x,r}:x\in\mathcal{X},r\geq
0\\}$. The VC dimension of $\mathcal{F}_{\mathcal{X}}$ is upper bounded by the
VC dimension of the set of all closed Euclidean balls in $\mathbb{R}^{d}$,
which is $d+1$; see [Dudley, 1979]. Therefore,

$\mathbb{E}\left[\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})\right]\leq\mathbb{E}\left[\sup_{\phi\in\mathcal{F}_{\mathcal{X}}}|\mu\phi-\hat{\mu}\phi|\right]\leq
C\sqrt{\frac{\mathrm{vc}(\mathcal{F}_{\mathcal{X}})}{n}}=C\sqrt{\frac{d+1}{n}}.$
(A.6)

Repeating the argument in the proof of Lemma 2, we can deduce that

$\sup_{\phi\in\mathcal{F}_{\mathcal{X}}}|\mu\phi-\hat{\mu}\phi|\leq\mathbb{E}\left[\sup_{\phi\in\mathcal{F}_{\mathcal{X}}}|\mu\phi-\hat{\mu}\phi|\right]+\sqrt{\frac{\log(1/\delta)}{2n}}$
(A.7)

holds with probability at least $1-\delta$. We obtain (3.1) by combining (A.7)
with (A.5) and (A.6). ∎

Next, we prove Theorem 2.

###### Proof of Theorem 2.

As the supports of $\mu_{x}$ and $\hat{\mu}_{x}$ are contained in the interval
$[0,\Delta]$ for any $x\in\mathcal{X}$, by the Kantorovich-Rubinstein theorem
(Theorem 1.14 of [Villani, 2003]), we have

$W_{1}(\mu_{x},\hat{\mu}_{x})=\sup_{f\in\mathrm{Lip}_{1}}\left(\mu_{x}f-\hat{\mu}_{x}f\right),$

where
$\mathrm{Lip}_{1}=\left\\{f\colon[0,\Delta]\to\mathbb{R}\,|\,f(0)=0~{}\text{and}~{}f~{}\text{is}~{}1\text{-Lipschitz}\right\\}$.
Hence,

$\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})=\sup_{x\in\mathcal{X}}\sup_{f\in\mathrm{Lip}_{1}}\left(\mu_{x}f-\hat{\mu}_{x}f\right)=\sup_{\phi\in\mathcal{F}}\left(\mu\phi-\hat{\mu}\phi\right),$

where $\mathcal{F}=\\{f\circ
d_{\mathcal{X},x}:f\in\mathrm{Lip}_{1},x\in\mathcal{X}\\}$ which is a function
class from $\mathcal{X}$ to $\mathbb{R}$; here, we are using
$d_{\mathcal{X},x}(\cdot)=d_{\mathcal{X}}(x,\cdot)$ as in Definition 5. Using
symmetrization and Dudley’s chaining argument (see Chapter 8 of [Vershynin,
2018]), one can derive the following:

$\mathbb{E}\left[\sup_{x\in\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x})\right]=\mathbb{E}\left[\sup_{\phi\in\mathcal{F}}\left(\mu\phi-\hat{\mu}\phi\right)\right]\leq\frac{C}{\sqrt{n}}\int_{0}^{\Delta/2}\sqrt{\log\mathcal{N}(\varepsilon,\mathcal{F},\|\cdot\|_{\mathcal{X}})}\,\mathrm{d}\varepsilon,$
(A.8)

where $\mathcal{N}(\varepsilon,\mathcal{F},\|\cdot\|_{\mathcal{X}})$ denotes
the $\varepsilon$-covering number of $\mathcal{F}$ with respect to
$\|\cdot\|_{\mathcal{X}}$, namely, the sup norm defined by
$\|\phi\|_{\mathcal{X}}=\sup_{x\in\mathcal{X}}|\phi(x)|$ for any
$\phi\colon\mathcal{X}\to\mathbb{R}$. Here, we are using the fact that
$\sup_{\phi\in\mathcal{F}}\|\phi\|_{\mathcal{X}}\leq\Delta$ since
$\sup_{z\in[0,\Delta]}|f(z)|\leq\Delta$ for any $f\in\mathrm{Lip}_{1}$ by
Lipschitzness. Now, noticing that $\mathcal{F}$ is a function class consisting
of composite functions, we upper bound
$\mathcal{N}(\varepsilon,\mathcal{F},\|\cdot\|_{\mathcal{X}})$ by the product
of the $\varepsilon$-covering number of $\mathrm{Lip}_{1}$ with respect to the
sup norm $\|\cdot\|_{\infty}$ on $[0,\Delta]$, that is,
$\|f\|_{\infty}=\sup_{z\in[0,\Delta]}|f(z)|$, and the $\varepsilon$-covering
number of $\mathcal{X}$ with respect to $d_{\mathcal{X}}$. More precisely, we
claim

$\mathcal{N}(2\varepsilon,\mathcal{F},\|\cdot\|_{\mathcal{X}})\leq\mathcal{N}(\varepsilon,\mathrm{Lip}_{1},\|\cdot\|_{\infty})\times\mathcal{N}(\varepsilon,\mathcal{X},d_{\mathcal{X}}).$
(A.9)

To see this, let $\mathrm{Lip}_{1}^{\varepsilon}\subset\mathrm{Lip}_{1}$ and
$\mathcal{X}^{\varepsilon}\subset\mathcal{X}$ be $\varepsilon$-coverings of
$\mathrm{Lip}_{1}$ and $\mathcal{X}$, respectively. For $f\in\mathrm{Lip}_{1}$
and $x\in\mathcal{X}$, we can find
$f^{\prime}\in\mathrm{Lip}_{1}^{\varepsilon}$ and
$x^{\prime}\in\mathcal{X}^{\varepsilon}$ such that
$\|f-f^{\prime}\|_{\infty}\leq\varepsilon$ and
$d_{\mathcal{X}}(x,x^{\prime})\leq\varepsilon$. Then,

$\begin{split}\|f\circ d_{\mathcal{X},x}-f^{\prime}\circ
d_{\mathcal{X},x^{\prime}}\|_{\mathcal{X}}&\leq\|f\circ
d_{\mathcal{X},x}-f\circ d_{\mathcal{X},x^{\prime}}\|_{\mathcal{X}}+\|f\circ
d_{\mathcal{X},x^{\prime}}-f^{\prime}\circ
d_{\mathcal{X},x^{\prime}}\|_{\mathcal{X}}\\\
&\leq\|d_{\mathcal{X},x}-d_{\mathcal{X},x^{\prime}}\|_{\mathcal{X}}+\|f-f^{\prime}\|_{\infty}\\\
&\leq d_{\mathcal{X}}(x,x^{\prime})+\varepsilon\\\ &\leq
2\varepsilon.\end{split}$

Hence, we obtain (A.9). The standard upper bound on the covering number of
$\mathrm{Lip}_{1}$ (see Exercise 8.2.7 of [Vershynin, 2018]) tells

$\log\mathcal{N}(\varepsilon,\mathrm{Lip}_{1},\|\cdot\|_{\infty})\leq\frac{c\Delta}{\varepsilon},$
(A.10)

where $c$ is some absolute constant. Combining (A.8), (A.9), and (A.10), we
obtain (3.2). Lastly, using the fact
$\sup_{\phi\in\mathcal{F}}\|\phi\|_{\mathcal{X}}\leq\Delta$ again and applying
McDiarmid’s inequality, we have (3.3). ∎

#### Convergence under $\mathcal{D}_{1}$

Lastly, we provide details of Remark 6. First, let us show that the rate of
convergence of $\mathbb{E}\mathcal{D}_{1}(\mu,\hat{\mu})$ can be as slow as
that of $\mathbb{E}W_{1}(\mu,\hat{\mu})$. To see this, notice that

$\begin{split}\mathcal{D}_{1}(\mu,\hat{\mu})&=\inf_{\gamma\in\Pi(\mu,\hat{\mu})}\int_{\mathcal{X}\times\mathcal{X}}W_{1}(\mu_{x},\hat{\mu}_{x^{\prime}})\,\mathrm{d}\gamma(x,x^{\prime})\\\
&\leq\underbrace{\inf_{\gamma\in\Pi(\mu,\hat{\mu})}\int_{\mathcal{X}\times\mathcal{X}}W_{1}(\mu_{x},\mu_{x^{\prime}})\,\mathrm{d}\gamma(x,x^{\prime})}_{(*)}+\underbrace{\int_{\mathcal{X}}W_{1}(\mu_{x^{\prime}},\hat{\mu}_{x^{\prime}})\,\mathrm{d}\hat{\mu}(x^{\prime})}_{(**)},\end{split}$

where the inequality is due to the triangle inequality
$W_{1}(\mu_{x},\hat{\mu}_{x^{\prime}})\leq
W_{1}(\mu_{x},\mu_{x^{\prime}})+W_{1}(\mu_{x^{\prime}},\hat{\mu}_{x^{\prime}})$.
To upper bound the expectation of (**), note that

$\mathbb{E}(**)=\mathbb{E}\int_{\mathcal{X}}W_{1}(\mu_{x^{\prime}},\hat{\mu}_{x^{\prime}})\,\mathrm{d}\hat{\mu}(x^{\prime})=\mathbb{E}\left(\frac{1}{n}\sum_{i=1}^{n}W_{1}(\mu_{X_{i}},\hat{\mu}_{X_{i}})\right)=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}W_{1}(\mu_{X_{i}},\hat{\mu}_{X_{i}}),$

which can be shown to converge at the rate of $O(n^{-1/2})$ under mild
assumptions by Lemma 3 below. As $W_{1}(\mu_{x},\mu_{x^{\prime}})\leq
d_{\mathcal{X}}(x,x^{\prime})$ by Lemma 1, we can see that (*) $\leq
W_{1}(\mu,\hat{\mu})$, meaning that the rate of convergence of
$\mathbb{E}\mathcal{D}_{1}(\mu,\hat{\mu})$ can be as slow as that of
$\mathbb{E}W_{1}(\mu,\hat{\mu})$. In the presence of symmetry, however, (*)
can be much smaller than $W_{1}(\mu,\hat{\mu})$. Importantly, consider the
rotational invariance, namely, suppose $\mathcal{X}=\mathbb{R}^{d}$ and $\mu$
is rotationally invariant, then $\mu_{x}=\mu_{x^{\prime}}$ if
$\|x\|_{2}=\|x^{\prime}\|_{2}$. Accordingly, for any
$x,x^{\prime}\in\mathbb{R}^{d}$,

$W_{1}(\mu_{x},\mu_{x^{\prime}})=W_{1}(\mu_{\|x\|_{2}e_{1}},\mu_{\|x^{\prime}\|_{2}e_{1}})\leq\left\|\|x\|_{2}e_{1}-\|x^{\prime}\|_{2}e_{1}\right\|_{2}=\left|\|x\|_{2}-\|x^{\prime}\|_{2}\right|,$
(A.11)

where $e_{1}=(1,0,\ldots,0)\in\mathbb{R}^{d}$. Now, define
$\Phi\colon\mathbb{R}^{d}\to\mathbb{R}$ by $\Phi(x)=\|x\|_{2}$. For any
coupling $\gamma\in\Pi(\mu,\hat{\mu})$, note that
$(\Phi\times\Phi)_{\\#}\gamma\in\Pi(\Phi_{\\#}\mu,\Phi_{\\#}\hat{\mu})$, where
$\Phi\times\Phi\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}\times\mathbb{R}$
is a map defined by
$(\Phi\times\Phi)(x,x^{\prime})=(\Phi(x),\Phi(x^{\prime}))$. Conversely, for
any coupling $\rho\in\Pi(\Phi_{\\#}\mu,\Phi_{\\#}\hat{\mu})$, we can find a
coupling $\gamma\in\Pi(\mu,\hat{\mu})$ such that
$\rho=(\Phi\times\Phi)_{\\#}\gamma$ by Lemma 3.12 of [Ambrosio and Gigli,
2013]. In summary,

$\Pi(\Phi_{\\#}\mu,\Phi_{\\#}\hat{\mu})=\\{(\Phi\times\Phi)_{\\#}\gamma:\gamma\in\Pi(\mu,\hat{\mu})\\}.$
(A.12)

Therefore,

$\begin{split}\inf_{\gamma\in\Pi(\mu,\hat{\mu})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}W_{1}(\mu_{x},\mu_{x^{\prime}})\,\mathrm{d}\gamma(x,x^{\prime})&=\inf_{\gamma\in\Pi(\mu,\hat{\mu})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}W_{1}(\mu_{\Phi(x)e_{1}},\mu_{\Phi(x^{\prime})e_{1}})\,\mathrm{d}\gamma(x,x^{\prime})\\\
&=\inf_{\gamma\in\Pi(\mu,\hat{\mu})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}W_{1}(\mu_{te_{1}},\mu_{se_{1}})\,\mathrm{d}(\Phi\times\Phi)_{\\#}\gamma(t,s)\\\
&=\inf_{\rho\in\Pi(\Phi_{\\#}\mu,\Phi_{\\#}\hat{\mu})}\int_{\mathbb{R}\times\mathbb{R}}W_{1}(\mu_{te_{1}},\mu_{se_{1}})\,\mathrm{d}\rho(t,s)\\\
&\leq\inf_{\rho\in\Pi(\Phi_{\\#}\mu,\Phi_{\\#}\hat{\mu})}\int_{\mathbb{R}\times\mathbb{R}}|t-s|\,\mathrm{d}\rho(t,s)\\\
&=W_{1}(\Phi_{\\#}\mu,\Phi_{\\#}\hat{\mu}),\end{split}$

where the third equality is due to (A.12) and the inequality is due to (A.11).
As $\Phi_{\\#}\hat{\mu}$ is essentially the empirical measure of
$\Phi_{\\#}\mu$ which is defined on $\mathbb{R}$, we conclude that

$\mathbb{E}(*)\leq\mathbb{E}W_{1}(\Phi_{\\#}\mu,\Phi_{\\#}\hat{\mu})=O(n^{-1/2}).$

Therefore, in this case, we obtain that
$\mathbb{E}\mathcal{D}_{1}(\mu,\hat{\mu})$ is bounded by $O(n^{-1/2})$,
whereas $\mathbb{E}W_{1}(\mu,\hat{\mu})$ can be as slow as $O(n^{-1/d})$.

######  Lemma 3.

Let $(\mathcal{X},d_{\mathcal{X}},\mu)$ be a metric measure space and
$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ be the empirical measure
based on $X_{1},\ldots,X_{n}$ that are i.i.d. from $\mu$. Suppose
$\Delta=\sup_{x,x^{\prime}\in\mathcal{X}}d_{\mathcal{X}}(x,x^{\prime})<\infty$
and the support of $\mu$ is exactly $\mathcal{X}$. Then, for each
$i\in\\{1,\ldots,n\\}$,

$\mathbb{E}\left[W_{1}(\mu_{X_{i}},\hat{\mu}_{X_{i}})\right]\leq\frac{C}{\sqrt{n-1}}+\frac{\Delta+1}{n},$
(A.13)

where $C$ is the absolute constant in Lemma 2.

###### Proof of Lemma 3.

Without loss of generality, let $i=1$. Also, let
$\phi_{x,r}(\cdot)=1\\{d_{\mathcal{X}}(x,\cdot)\leq r\\}$ as in the proof of
Lemma 2. Then, for any $r\geq 0$,

$\hat{\mu}\phi_{X_{1},r}=\frac{1}{n}\sum_{k=1}^{n}\phi_{X_{1},r}(X_{k})=\frac{1}{n}\sum_{k=2}^{n}\phi_{X_{1},r}(X_{k})+\frac{1}{n}=\frac{n-1}{n}\hat{\mu}_{-1}\phi_{X_{1},r}+\frac{1}{n},$

where $\hat{\mu}_{-1}=\frac{1}{n-1}\sum_{k=2}^{n}\delta_{X_{k}}$. Therefore,

$\begin{split}|\mu\phi_{X_{1},r}-\hat{\mu}\phi_{X_{1},r}|&=\left|\mu\phi_{X_{1},r}-\hat{\mu}_{-1}\phi_{X_{1},r}+\frac{\hat{\mu}_{-1}\phi_{X_{1},r}}{n}-\frac{1}{n}\right|\\\
&\leq|\mu\phi_{X_{1},r}-\hat{\mu}_{-1}\phi_{X_{1},r}|+\frac{\hat{\mu}_{-1}\phi_{X_{1},r}}{n}+\frac{1}{n}.\end{split}$

Accordingly,

$\begin{split}W_{1}(\mu_{X_{1}},\hat{\mu}_{X_{1}})&=\int_{0}^{\Delta}|\mu\phi_{X_{1},r}-\hat{\mu}\phi_{X_{1},r}|\,\mathrm{d}r\\\
&\leq\int_{0}^{\Delta}|\mu\phi_{X_{1},r}-\hat{\mu}_{-1}\phi_{X_{1},r}|\,\mathrm{d}r+\frac{1}{n}\int_{0}^{\Delta}\hat{\mu}_{-1}\phi_{X_{1},r}\,\mathrm{d}r+\frac{1}{n}\\\
&\leq\sup_{\phi\in\mathcal{F}_{X_{1}}}|\mu\phi-\hat{\mu}_{-1}\phi|+\frac{1}{n}\int_{0}^{\Delta}\hat{\mu}_{-1}\phi_{X_{1},r}\,\mathrm{d}r+\frac{1}{n},\end{split}$

where $\mathcal{F}_{x}:=\\{\phi_{x,r}:r\geq 0\\}$ as in the proof of Lemma 2.
Since $\mathbb{E}[\hat{\mu}_{-1}\phi_{X_{1},r}|X_{1}]=\mu\phi_{X_{1},r}$, we
have

$\mathbb{E}\left[W_{1}(\mu_{X_{1}},\hat{\mu}_{X_{1}})|X_{1}\right]\leq\mathbb{E}\left[\sup_{\phi\in\mathcal{F}_{X_{1}}}|\mu\phi-\hat{\mu}_{-1}\phi|\,\Big{|}\,X_{1}\right]+\frac{1}{n}\int_{0}^{\Delta}\mu\phi_{X_{1},r}\,\mathrm{d}r+\frac{1}{n}.$

As in the proof of Lemma 2, we have

$\mathbb{E}\left[\sup_{\phi\in\mathcal{F}_{X_{1}}}|\mu\phi-\hat{\mu}_{-1}\phi|\,\Big{|}\,X_{1}\right]\leq
C\sqrt{\frac{\mathrm{vc}(\mathcal{F}_{X_{1}})}{n-1}}=\frac{C}{\sqrt{n-1}}.$

Meanwhile,

$\int_{0}^{\Delta}\mu\phi_{X_{1},r}\,\mathrm{d}r\leq\int_{0}^{\Delta}1\,\mathrm{d}r=\Delta.$

Hence,

$\mathbb{E}[W_{1}(\mu_{X_{1}},\hat{\mu}_{X_{1}})|X_{1}]\leq\frac{C}{\sqrt{n-1}}+\frac{\Delta}{n}+\frac{1}{n},$

which leads to (A.13) by the law of iterated expectations. ∎

###  A.2 Proof of Theorem 3

###### Proof of Theorem 3.

For the proof, we may assume $\pi^{\ast}$ is the identity without loss of
generality. First, letting
$Q(\pi)=\sum_{i=1}^{n}W_{1}(\hat{\mu}_{i},\hat{\nu}_{\pi(i)})$, observe that

$\left(\hat{\pi}\neq\pi^{\ast}\right)\subset\bigcup_{\pi\in\mathcal{S}_{n}}\left(Q(\pi^{\ast})\geq
Q(\pi)\right).$

Now, the event $Q(\pi^{\ast})\geq Q(\pi)$ is equivalent to

$\sum_{i=1}^{n}\left(W_{1}(\hat{\mu}_{i},\hat{\nu}_{i})-W_{1}(\hat{\mu}_{i},\hat{\nu}_{\pi(i)})\right)=\sum_{i:\pi(i)\neq
i}\left(W_{1}(\hat{\mu}_{i},\hat{\nu}_{i})-W_{1}(\hat{\mu}_{i},\hat{\nu}_{\pi(i)})\right)\geq
0,$

which is contained in the following events:

$\bigcup_{i:\pi(i)\neq i}\left(W_{1}(\hat{\mu}_{i},\hat{\nu}_{i})\geq
W_{1}(\hat{\mu}_{i},\hat{\nu}_{\pi(i)})\right)\subset\bigcup_{i=1}^{n}\bigcup_{j\neq
i}\left(W_{1}(\hat{\mu}_{i},\hat{\nu}_{i})\geq
W_{1}(\hat{\mu}_{i},\hat{\nu}_{j})\right)$

Therefore,

$\left(\hat{\pi}\neq\pi^{\ast}\right)\subset\bigcup_{i=1}^{n}\bigcup_{j\neq
i}\left(W_{1}(\hat{\mu}_{i},\hat{\nu}_{i})\geq
W_{1}(\hat{\mu}_{i},\hat{\nu}_{j})\right).$ (A.14)

Now, recall that we can write a rigid transformation
$T\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ as $T(\theta)=R\theta+b$ for some
orthonormal matrix $R\in\mathbb{R}^{d\times d}$ and some vector
$b\in\mathbb{R}^{d}$; accordingly,

$\displaystyle\hat{\mu}_{i}$
$\displaystyle:=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|X_{i}-X_{k}\|_{2}^{2}}=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|\theta_{i}-\theta_{k}+\sigma\xi_{ik}\|_{2}^{2}},$
$\displaystyle\hat{\nu}_{i}$
$\displaystyle:=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|Y_{i}-Y_{k}\|_{2}^{2}}=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|R(\theta_{i}-\theta_{k})+\sigma(\zeta_{i}-\zeta_{k})\|_{2}^{2}}=\frac{1}{n}\sum_{k=1}^{n}\delta_{\|\theta_{i}-\theta_{k}+\sigma\zeta_{ik}\|_{2}^{2}},$

where we let $\xi_{ik}=\xi_{i}-\xi_{k}$ and
$\zeta_{ik}=R^{\top}(\zeta_{i}-\zeta_{k})$. For $W_{1}$ between probability
measures supported on the same number of points in $\mathbb{R}$ with uniform
weights, we have the following:

$W_{1}(\hat{\mu}_{i},\hat{\nu}_{j})=\min_{\pi\in\mathcal{S}_{n}}\frac{1}{n}\sum_{k=1}^{n}\left|\|\theta_{i}-\theta_{k}+\sigma\xi_{ik}\|_{2}^{2}-\|\theta_{j}-\theta_{\pi(k)}+\sigma\zeta_{j\pi(k)}\|_{2}^{2}\right|.$

Note that,

$\begin{split}&\left|\|\theta_{i}-\theta_{k}+\sigma\xi_{ik}\|_{2}^{2}-\|\theta_{j}-\theta_{\pi(k)}+\sigma\zeta_{j\pi(k)}\|_{2}^{2}\right|\\\
&=\left|\|\theta_{i}-\theta_{k}\|_{2}^{2}-\|\theta_{j}-\theta_{\pi(k)}\|_{2}^{2}+\sigma^{2}(\|\xi_{ik}\|_{2}^{2}-\|\zeta_{j\pi(k)}\|_{2}^{2})+2\sigma\left((\theta_{i}-\theta_{k})^{\top}\xi_{ik}-(\theta_{j}-\theta_{\pi(k)})^{\top}\zeta_{j\pi(k)}\right)\right|\\\
&\geq\left|\|\theta_{i}-\theta_{k}\|_{2}^{2}-\|\theta_{j}-\theta_{\pi(k)}\|_{2}^{2}\right|-\sigma^{2}Z-2\sigma(U+V),\end{split}$

where we define

$\displaystyle Z$
$\displaystyle:=\max_{k,k^{\prime},\ell,\ell^{\prime}}\left|\|\xi_{k\ell}\|_{2}^{2}-\|\zeta_{k^{\prime}\ell^{\prime}}\|_{2}^{2}\right|,$
$\displaystyle U$
$\displaystyle:=\max_{k\neq\ell}|(\theta_{k}-\theta_{\ell})^{\top}\xi_{k\ell}|,$
$\displaystyle V$
$\displaystyle:=\max_{k\neq\ell}|(\theta_{k}-\theta_{\ell})^{\top}\zeta_{k\ell}|.$

Therefore,

$\begin{split}W_{1}(\hat{\mu}_{i},\hat{\nu}_{j})&\geq\min_{\pi\in\mathcal{S}_{n}}\frac{1}{n}\sum_{k=1}^{n}\left|\|\theta_{i}-\theta_{k}\|_{2}^{2}-\|\theta_{j}-\theta_{\pi(k)}\|_{2}^{2}\right|-\sigma^{2}Z-2\sigma(U+V)\\\
&=W_{1}(\mu_{i},\mu_{j})-\sigma^{2}Z-2\sigma(U+V)\\\
&\geq\Omega-\sigma^{2}Z-2\sigma(U+V).\end{split}$

Similarly, we can show $W_{1}(\hat{\mu}_{i},\hat{\nu}_{j})\leq
W_{1}(\mu_{i},\mu_{j})+\sigma^{2}Z+2\sigma(U+V)$. Accordingly, we have

$\left|W_{1}(\hat{\mu}_{i},\hat{\nu}_{j})-W_{1}(\mu_{i},\mu_{j})\right|\leq\sigma^{2}Z+2\sigma(U+V)\quad\forall
i,j\in\\{1,\ldots,n\\}.$

In particular, for any $i\in\\{1,\ldots,n\\}$, we have

$W_{1}(\hat{\mu}_{i},\hat{\nu}_{i})\leq\sigma^{2}Z+2\sigma(U+V).$

Therefore, if we can find $i,j\in\\{1,\ldots,n\\}$ such that $i\neq j$ and
$W_{1}(\hat{\mu}_{i},\hat{\nu}_{i})\geq W_{1}(\hat{\mu}_{i},\hat{\nu}_{j})$,
we must have

$\sigma^{2}Z+2\sigma(U+V)\geq\Omega-\sigma^{2}Z-2\sigma(U+V).$

Accordingly, by (A.14), we have

$\mathbb{P}(\hat{\pi}\neq\pi^{\ast})\leq\mathbb{P}(2\sigma^{2}Z+4\sigma(U+V)\geq\Omega)$

Now, let

$Z=\max_{k,k^{\prime},\ell,\ell^{\prime}}\left|\|\xi_{k\ell}\|_{2}^{2}-\|\zeta_{k^{\prime}\ell^{\prime}}\|_{2}^{2}\right|\leq
2\max_{k,\ell}\left|\frac{\|\xi_{k\ell}\|_{2}^{2}}{2}-d\right|+2\max_{k,\ell}\left|\frac{\|\zeta_{k\ell}\|_{2}^{2}}{2}-d\right|=:2Z_{1}+2Z_{2}.$

Then,

$\mathbb{P}(\hat{\pi}\neq\pi^{\ast})\leq\mathbb{P}\left(Z_{1}\geq\frac{\Omega}{16\sigma^{2}}\right)+\mathbb{P}\left(Z_{2}\geq\frac{\Omega}{16\sigma^{2}}\right)+\mathbb{P}\left(U\geq\frac{\Omega}{16\sigma}\right)+\mathbb{P}\left(V\geq\frac{\Omega}{16\sigma}\right).$

We use the following concentration inequality of the chi-squared random
variable from Lemma 10 of [Collier and Dalalyan, 2016]:

$\mathbb{P}(|\chi^{2}(d)-d|\geq t)\leq
2\exp\left(-\frac{t}{8}\min\left(\frac{t}{d},1\right)\right)\quad\forall t>0.$

As $\|\xi_{k\ell}\|_{2}^{2},\|\zeta_{k\ell}\|^{2}\sim 2\chi^{2}(d)$ for
$k\neq\ell$, we have

$\begin{split}\mathbb{P}\left(Z_{1}\geq\frac{\Omega}{16\sigma^{2}}\right)&=\mathbb{P}\left(Z_{2}\geq\frac{\Omega}{16\sigma^{2}}\right)\\\
&\leq\sum_{k\neq\ell}\mathbb{P}\left(\left|\frac{\|\zeta_{k\ell}\|_{2}^{2}}{2}-d\right|\geq\frac{\Omega}{16\sigma^{2}}\right)\\\
&\leq
2n^{2}\exp\left(-\frac{\Omega}{128\sigma^{2}}\min\left(\frac{\Omega}{16\sigma^{2}d},1\right)\right).\end{split}$

Meanwhile, as
$(\theta_{k}-\theta_{\ell})^{\top}\xi_{k\ell},(\theta_{k}-\theta_{\ell})^{\top}\zeta_{k\ell}\sim
N(0,2\|\theta_{k}-\theta_{\ell}\|_{2}^{2})$,

$\begin{split}\mathbb{P}\left(U\geq\frac{\Omega}{16\sigma}\right)&=\mathbb{P}\left(V\geq\frac{\Omega}{16\sigma}\right)\\\
&\leq\sum_{k\neq\ell}\mathbb{P}\left(|(\theta_{k}-\theta_{\ell})^{\top}\zeta_{k\ell}|\geq\frac{\Omega}{16\sigma}\right)\\\
&\leq\sum_{k\neq\ell}2\exp\left(-\frac{\Omega^{2}}{1024\sigma^{2}\|\theta_{k}-\theta_{\ell}\|_{2}^{2}}\right)\\\
&\leq
2n^{2}\exp\left(-\frac{\Omega^{2}}{1024\sigma^{2}\max_{k,\ell}\|\theta_{k}-\theta_{\ell}\|_{2}^{2}}\right).\end{split}$

Therefore,

$\begin{split}\mathbb{P}(\hat{\pi}\neq\pi^{\ast})&\leq
4n^{2}\exp\left(-\frac{\Omega}{128\sigma^{2}}\min\left(\frac{\Omega}{16\sigma^{2}d},1\right)\right)+4n^{2}\exp\left(-\frac{\Omega^{2}}{1024\sigma^{2}\Theta^{2}}\right)\\\
&\leq
8n^{2}\exp\left(-\min\left(\frac{\Omega}{128\sigma^{2}}\min\left(\frac{\Omega}{16\sigma^{2}d},1\right),\frac{\Omega^{2}}{1024\sigma^{2}\Theta^{2}}\right)\right)\\\
&=8n^{2}\exp\left(-\frac{\Omega}{128\sigma^{2}}\min\left(\frac{\Omega}{16\sigma^{2}d},\frac{\Omega}{8\Theta^{2}},1\right)\right).\end{split}$

Hence, in order to have $\mathbb{P}(\hat{\pi}\neq\pi^{\ast})\leq\delta$, we
need

$\log(8n^{2}/\delta)\leq\frac{\Omega}{128\sigma^{2}}\min\left(\frac{\Omega}{16\sigma^{2}d},\frac{\Omega}{8\Theta^{2}},1\right).$

Note that $W_{1}(\mu_{i},\mu_{j})\leq\Theta^{2}$, implying
$\Omega\leq\Theta^{2}$. Hence, it suffices to have

$\log(8n^{2}/\delta)\leq\frac{\Omega}{128\sigma^{2}}\min\left(\frac{\Omega}{16\sigma^{2}d},\frac{\Omega}{8\Theta^{2}}\right).$

Equivalently,

$1024\sigma^{2}\log(8n^{2}/\delta)\cdot\max\left(\Theta^{2},2\sigma^{2}d\right)\leq\Omega^{2}.$

∎