Fr\'echet Means and Procrustes Analysis in Wasserstein Space

TL;DR

This paper explores the computation of Fréchet means and optimal registration of measures in Wasserstein space, linking these problems through a geometric perspective and proposing consistent estimators for practical data analysis.

Contribution

It introduces a novel connection between Fréchet means and Procrustes analysis in Wasserstein space, with methods for their computation and statistical inference.

Findings

Proposed a gradient descent method for Fréchet mean computation.

Established the equivalence between Fréchet mean and Procrustes registration.

Proved consistency of nonparametric estimators for population mean and registration maps.

Abstract

We consider two statistical problems at the intersection of functional and non-Euclidean data analysis: the determination of a Fr\'echet mean in the Wasserstein space of multivariate distributions; and the optimal registration of deformed random measures and point processes. We elucidate how the two problems are linked, each being in a sense dual to the other. We first study the finite sample version of the problem in the continuum. Exploiting the tangent bundle structure of Wasserstein space, we deduce the Fr\'echet mean via gradient descent. We show that this is equivalent to a Procrustes analysis for the registration maps, thus only requiring successive solutions to pairwise optimal coupling problems. We then study the population version of the problem, focussing on inference and stability: in practice, the data are i.i.d. realisations from a law on Wasserstein space, and indeed their observation is discrete, where one observes a proxy finite sample or point process. We construct regularised nonparametric estimators, and prove their consistency for the population mean, and uniform consistency for the population Procrustes registration maps.

Figures

Equations

$${\int }_{{\mathbb{R}}^d}{\Vert {\mathbf{t}}_{\mu }^{\nu }(x)-x\Vert }^2\mathrm{d}\mu (x)\le {\int }_{{\mathbb{R}}^d}{\Vert \mathbf{q}(x)-x\Vert }^2\mathrm{d}\mu (x),$$

$${\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\Vert x-y{\Vert }^2\mathrm{d}\xi (x,y)$$

$$\mathbb{E}\Vert X-Y{\Vert }^2,$$

$$d(\mu ,\nu )=\underset{\xi \in \Gamma (\mu ,\nu )}{\inf }\sqrt{{\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\Vert x-y{\Vert }^2\mathrm{d}\xi (x,y)},$$

$$d(\mu ,{\mu }^i)=\sqrt{{\int }_{{\mathbb{R}}^d}{\Vert {\mathbf{t}}_{\mu }^{{\mu }^i}(x)-x\Vert }^2\mathrm{d}\mu (x)}.$$

$${\mathbf{t}}_{\mu }^{\nu }=G_{\nu }^{-1}\circ G_{\mu },$$

$$\mathbf{t}=\nabla \phi .$$

$$F(\gamma ):=\frac{1}{2N}\sum _{i=1}^N{d}^2({\mu }^i,\gamma )$$

$$\mathbf{X}=(X_1,\dots ,X_N),X_i={\mathbf{t}}_{\bar{\mu }}^{{\mu }^i}(Z),i=1\dots ,N,$$

$$\sum _{i=1}^N\sum _{j=i+1}^N\mathbb{E}\Vert X_i-X_j{\Vert }^2\le \sum _{i=1}^N\sum _{j=i+1}^N\mathbb{E}\Vert Y_i-Y_j{\Vert }^2$$

$${\mathbf{t}}_{\bar{\mu }}^{{\mu }^i}\#\bar{\mu }={\mu }^i.$$

$${\mu }_t=\left[\mathbf{i}+t({\mathbf{t}}_{{\mu }^0}^{{\mu }^1}-\mathbf{i})\right]\#{\mu }^0,t\in [0,1].$$

$${\mathrm{Tan}}_{\mu }={\mathrm{Tan}}_{\mu }{\mathcal{P}}_2({\mathbb{R}}^d)={\overline{\{\nabla \phi :\phi \in C_c^{\infty }({\mathbb{R}}^d)\}}}^{L^2(\mu )},$$

$${\mathrm{Tan}}_{\mu }={\overline{\{\lambda (\mathbf{r}-\mathbf{i}):\mathbf{r}\text{ optimal between }\mu \text{ and }\mathbf{r}\#\mu ;\lambda >0\}}}^{L^2(\mu )},$$

$${\exp }_{{\mu }^0}(\mathbf{r}-\mathbf{i})=\mathbf{r}\#{\mu }^0$$

$${\log }_{{\mu }^0}({\mu }^1)={\mathbf{t}}_{{\mu }^0}^{{\mu }^1}-\mathbf{i},$$

$$F_0:{\mathcal{P}}_2({\mathbb{R}}^d)\to \mathbb{R},F_0(\mu )=\frac{1}{2}{d}^2(\mu ,{\mu }^0).$$

$$\underset{\nu \to \mu }{\lim }\frac{F_0(\nu )-F_0(\mu )+{\int }_{{\mathbb{R}}^d}⟨{\mathbf{t}}_{\mu }^{{\mu }^0}(x)-x,{\mathbf{t}}_{\mu }^{\nu }(x)-x⟩\mathrm{d}\mu (x)}{d(\nu ,\mu )}=0,$$

$$⟨{\mathbf{t}}_{\mu }^{{\mu }^0}-\mathbf{i},{\mathbf{t}}_{\mu }^{\nu }-\mathbf{i}⟩$$

$$F_0(\nu )-F_0(\mu )=-⟨{\log }_{\mu }({\mu }^0),{\log }_{\mu }(\nu )⟩+o(d(\nu ,\mu )),\nu \to \mu ,$$

$$F_0^{\prime }(\mu )=-{\log }_{\mu }({\mu }^0)=-\left({\mathbf{t}}_{\mu }^{{\mu }^0}-\mathbf{i}\right)\in {\mathrm{Tan}}_{\mu }.$$

$$F(\gamma )=\frac{1}{2N}\sum _{i=1}^N{d}^2(\gamma ,{\mu }^i),\gamma \in {\mathcal{P}}_2({\mathbb{R}}^d).$$

$$F^{\prime }(\gamma )=-\frac{1}{N}\sum _{i=1}^N{\log }_{\gamma }({\mu }^i)=-\frac{1}{N}\sum _{i=1}^N\left({\mathbf{t}}_{\gamma }^{{\mu }^i}-\mathbf{i}\right).$$

$$\frac{1}{N}\sum _{i=1}^N\left({\mathbf{t}}_{\mu }^{{\mu }^i}-\mathbf{i}\right)=0,\mu -\text{mbox}almosteverywhere.$$

$${\gamma }_{j+1}={\exp }_{{\gamma }_j}\left(-{\tau }_j{F}^{\prime }({\gamma }_j)\right)=\left[\mathbf{i}+{\tau }_j\frac{1}{N}\sum _{i=1}^N{\log }_{\gamma }({\mu }^i)\right]\#{\gamma }_j=\left[\mathbf{i}+{\tau }_j\frac{1}{N}\sum _{i=1}^N({\mathbf{t}}_{{\gamma }_j}^{{\mu }^i}-\mathbf{i})\right]\#{\gamma }_j.$$

$$F({\gamma }_1)-F({\gamma }_0)\le -\Vert F^{\prime }({\gamma }_0){\Vert }^2\left[\tau -\frac{{\tau }^2}{2}\right].$$

$$d({\gamma }_j,\bar{\mu })\stackrel{j\to \infty }{⟶}0.$$

$$\underset{{\Omega }_1}{\sup }\Vert {\mathbf{t}}_{{\gamma }_j}^{{\mu }^i}-{\mathbf{t}}_{\bar{\mu }}^{{\mu }^i}\Vert \stackrel{j\to \infty }{⟶}0,\underset{{\Omega }_2}{\sup }\Vert {\mathbf{t}}_{{\mu }^i}^{{\gamma }_j}-{\mathbf{t}}_{{\mu }^i}^{\bar{\mu }}\Vert \stackrel{j\to \infty }{⟶}0,i=1,\dots ,N,$$

$$\left({\mathbf{t}}_{{\gamma }_j}^{{\mu }^1},\dots {\mathbf{t}}_{{\gamma }_j}^{{\mu }^n}\right)\#{\gamma }_j$$

$$\mathbb{E}[{\Pi }_i(A)]=\lambda (A),$$

Full text

62M30; 60D05,

62G07; 60G55,

Functional Data Analysis,

Manifold Statistics,

Optimal Transportation,

Phase Variation,

Point Process,

Random Measure,

Registration,

Warping,

and

T1Research supported by an ERC Starting Grant Award to Victor M. Panaretos.

1. 1 Introduction

2. 2 Optimal Transportation and Wasserstein Space

3. 3 Sample Setting

4. 3.1 Fréchet Means and Optimal Registration

5. 3.2 Wasserstein Geometry and the Gradient of the Fréchet Functional

6. 3.3 Gradient Descent and Procrustes Analysis

7. 4 Population Setting

8. 4.1 Deformation Models and Discrete Observation

9. 4.2 Regularised Nonparametric Estimation

10. 4.3 Asymptotic Theory

11. 5 Proofs of Formal Statements

12. 5.1 Proofs of Statements in Section 3.1

13. 5.2 Proofs of Statements in Section 3.2

14. 5.3 Proofs of Statements in Section 3.3

15. 5.4 Proofs of Statements in Section 4.1

16. 5.5 Monotone Operators, Optimal Transportation, Stochastic Convergence

17. 6 Some Examples

18. 6.1 The case $d=1$

19. 6.2 Independence

20. 6.3 Common Copulas

21. 6.4 Gaussian measures

22. 6.5 Partially Gaussian Trivariate Measures

23. 7 Concluding Remarks

24. 8 Supplementary Material

1 Introduction

Functional data analysis (e.g. Hsing & Eubank [41]) and non-Euclidean statistics (e.g. Patrangenaru & Ellingson [63]) represent modern areas of statistical research, whose key challenges arise from the intrinsic complexity of the data and the peculiarities of their ambient space. In the first case, the data are random elements in a separable Hilbert space of functions (typically $L^{2}[0,1]$), and resulting challenges are linked to infinite dimensionality (e.g. ill-posed studentisation, Munk et al. [57], and discrete measurements of continuum random objects, Zhang & Wang [78]). In the second case, the data are seen as random elements of a finite-dimensional Riemannian manifold (often a shape space), and resulting challenges are linked to the non-linear structure of the space (e.g. existence/uniqueness of Fréchet means, Le [51] and Kendall [46], and analysis of manifold variation, Huckemann, Munk & Hotz [42]).

At the intersection of these two domains, with manifestations in neurophysiology, imaging, and environmetrics, one finds data objects that are best modelled as distributions over $\mathbb{R}^{d}$, that is, random measures (Stoyan, Kendall & Mecke [24], Kallenberg [44]). Such random measures carry the infinite dimensional traits of functional data, but at the same time are characterised by intrinsic non-linearities due to their positivity and integrability constraints, requiring a non-Euclidean point of view. Indeed, despite their functional nature, their dominating variational feature is not due to additive amplitude fluctuations (as can be seen in the Karhunen-Loève expansion of functional data), but rather to random deformation of a structural mean (as in Freitag & Munk [34]) or template (as in morphometrics, Bookstein [20]). Still, being infinite dimensional, their observation is typically done discretely, for example noisily over a grid (e.g. Amit et al. [8], Allassonnière et al. [4]) or via random sampling (e.g. Panaretos & Zemel [61]), requiring tools and techniques from nonparametric statistics, as used in functional data analysis.

In this setting, the typical statistical objective is to estimate the underlying template that gives rise to the data by random deformation. This can often be modelled as a Fréchet mean with respect to some metric structure; dual to this problem is the recovery the deformation maps themselves, in order to register the individual realisations in a common coordinate system, given by registration maps. These problems are interwoven in shape theory, where the template and registration maps are the two ingredients of Procrustes analysis (Gower [38]; Dryden & Mardia [30]) and non-Euclidean PCA (Huckemann, Munk & Hotz [42]; Huckemann & Ziezold [43]). Obviously, the methods and algorithms for estimating a mean and carrying out a registration/Procrustes analysis are inextricably linked with the geometry of the sample space, which can be a matter of modelling choice or of first principles.

In this paper, we choose to study the problem of Fréchet averaging and Procrustes registration when the data are viewed as elements of the $L^{2}$-Wasserstein space of multivariate measures on $\mathbb{R}^{d}$. We choose this setting since it has a long history in assessing compatibility and fit of distributions related via deformations (Munk & Czado [56]; Freitag & Munk [34]), and as it can be seen to be a natural analogue of using $L^{2}$, in the case of measures111In the sense that the Wasserstein space is topologically homeomorphic to a convex subset of $L^{2}([0,1]^{d})$; when $d=1$, this homeomorphism is an isometry, whereas for $d>1$, it is a local isometry. (Panaretos & Zemel [61]; Bigot & Klein [14]). We work at both a sample level and a population level, as well as both at the level of continuum and discrete observation: our object of study is the determination of the Fréchet mean and registration maps at the level of a sample, as well as at their estimation when the observed measures are discretely observed realisations from a population of random measures. When $d=1$, the problem is well understood, owing to the flat geometry of Wasserstein space (Panaretos & Zemel [61]). When $d>1$, however, the Wasserstein space has non-negative curvature, and one encounters the classical difficulties of non-Euclidean statistics, augmented by the infinite dimensionality and discrete measurement of the problem (see Anderes et al. [9], Sommerfeld & Munk [71] and Tameling et al. [73] for challenges involved in the discrete setting).

In more detail, our contributions are:

• (A)

At the sample level: we illustrate how knowledge of the Fréchet mean (template) gives an explicit solution to the optimal registration/multicoupling problem (Section 3.1, Proposition 2). We study the tangent space geometry, using it to determine the gradient of the Fréchet functional (Section 3.2.2, Theorem 1), and characterise Karcher means via its zeroes (Corollary 1, Section 3.2.3). We give criteria for determining when a Karcher mean (local optimum) is a Fréchet mean (global optimum; Theorem 2). We construct a gradient descent algorithm (Algorithm 1), and find its optimal stepsize (Lemma 2) illustrating the algorithm structurally equivalent to a Procrustes algorithm (Section 3.3), reducing the determination of the mean to the successive solution of pairwise optimal transport problems. We prove that the gradient iterate converges to a Karcher mean in the Wasserstein metric (Subsection 3.3.2, Theorem 3); and that the induced transportation maps converge uniformly to the Procrustes maps (required for optimal mutlicoupling; Theorem 4, Section 3.3.3). The latter is particularly involved and requires techniques from the geometry of monotone operators on $\mathbb{R}^{d}$. As a noteworthy corollary, we deduce convergence of the multicouplings (Corollary 3).

• (B)

At the population level: we consider a population level model linking Fréchet means and optimal registration and give conditions for model identifiability (Section 4.1, Theorem 5); We then tackle the problem of point estimation of the population mean and registration maps in a functional data analysis setup, where instead of observing an i.i.d. sample $\{\mu^{1},\dots,\mu^{N}\}$ from the population, we observe samples or point processes with these measures as distributions/intensities. In this setting, we construct regularised nonparametric estimators of the Fréchet means and Procrustes maps, and prove that they are consistent in Wasserstein distance and uniform norm, respectively (Theorems 6 and 7).

Before presenting our main results, we first provide a short introduction to Wasserstein space in Section 2. Section 5 gathers the main proofs, for the sake of tidiness, and Section 6 presents several interesting special examples as an illustration. Section 8 supplements the main article, providing further technical details.

In reviewing an earlier version of our paper ([77], February 2016), a referee brought to our attention independent parallel work by Álvarez-Esteban et al., that had concurrently (January 2016) been submitted for publication in an analysis journal (and has now appeared, see [6]). Their work overlaps with part of ours in (A) above (Subsections 3.3.1 and 3.3.2). In particular, they too arrive at a (structurally) same algorithm (Algorithm 1). Their motivation, construction, and convergence proof differ substantially from ours (theirs is a fixed point iteration heuristically motivated by the Gaussian case, while their proof uses almost sure representations). Indeed, our geometrical framework and proof techniques is what allows us to study the problem of optimal registration (Procrustes analysis), requiring a careful study of the stochastic convergence of monotone operators on $\mathbb{R}^{d}$ (Section 5.5).

2 Optimal Transportation and Wasserstein Space

The reason the Wasserstein space arises as the natural space to capture deformation-based variation of random measures lies in its deep connection with the problem of optimal transportation of measure. This consists in solving the Monge problem (Villani [74]): given a pair of measures $(\mu,\nu)$, find a mapping $\mathbf{t}_{\mu}^{\nu}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ such that $\mathbf{t}_{\mu}^{\nu}\#\mu=\nu$, and

[TABLE]

for any other $\mathbf{q}$ such that $\mathbf{q}\#\mu=\nu$. Here, “$\#$” denotes the push-forward operation, where $[\mathbf{t}\#\mu](A)=\mu(\mathbf{t}^{-1}(A))$ for all Borel sets $A$ of $\mathbb{R}^{d}$. The map $\mathbf{t}_{\mu}^{\nu}$ is called an optimal transport plan, and a solution to this problem yields an optimal deformation of $\mu$ into $\nu$ with respect to the transport cost given by squared Euclidean distance.

An optimal transport map may fail to exist, and instead, one may need to solve the relaxed Monge problem, known as the Kantorovich problem (Villani [74]). Here instead of seeking a map $\mathbf{t}_{\mu}^{\nu}\#\mu=\nu$, one seeks a distribution $\xi$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $\mu$ and $\nu$, minimising the functional

[TABLE]

over all measures $\xi$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $\mu$ and $\nu$. In probabilistic terms, $\xi$ yields a coupling of random variables $X\sim\mu$ and $Y\sim\nu$ that minimises the quantity

[TABLE]

over all possible couplings of $X$ and $Y$. It can be shown that when the measure $\mu$ is regular (absolutely continuous with respect to Lebesgue measure), the Kantorovich problem reduces to the Monge problem, and the optimal coupling $\xi$ is supported on the graph of the function. That is, the optimal coupling exists, is unique, and can be realised by a proper transport map $\mathbf{t}_{\mu}^{\nu}$.

One may consider the space $\mathcal{P}_{2}(\mathbb{R}^{2})$ of all probability measures $\mu$ on $\mathbb{R}^{d}$ with finite variance (that is, ${\int_{\mathbb{R}^{d}}\!\|x\|^{2}\,\mathrm{d}\mu(x)}<\infty$) as a metric space, endowed with the $L^{2}$-Wasserstein distance

[TABLE]

where $\Gamma(\mu,\nu)$ is the set of probability measures on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $\mu$ and $\nu$. The induced metric space is colloquially called Wasserstein space and will form the geometrical context for our study of deformation-based variation of random measures. This space has been used extensively in statistics, as it metrises the topology of weak convergence, and convergence with respect to the metric yields both convergence in law, as well as convergence of the first two moments (for instance, in applications to the bootstrap, e.g. Bickel & Freedman [12], and to goodness-of-fit, e.g. Rippl, Munk & Sturm [65]).

The appropriateness of this distance when modeling deformations of measures becomes clear based on our previous remark concerning regularity: one can imagine an initial regular template $\mu$, that is deformed according to maps $\mathbf{q}_{i}$ to yield new measures $\mu^{i}=(\mathbf{q}_{i})\#\mu$. It is then natural to quantify the distance of the template to its perturbations by means of the minimal transportation (or deformation) cost

[TABLE]

That the distance can be expressed via a proper map, is due to the assumed regularity of $\mu$. Note that the maps $\mathbf{q}_{i}$ themselves will, in general, not be identifiable (many Borel maps can push $\mu$ forward to $\mu^{i}$). But they can be assumed to be exactly optimal, i.e. $\mathbf{q}_{i}=\mathbf{t}_{\mu}^{\mu^{i}}$ as a matter of parsimony, and in any case without loss of generality, leading to identifiability. These maps will also solve the registration problem: a map of the form $\mathbf{t}_{\mu}^{\mu^{i}}-\mathbf{i}$, with $\mathbf{i}$ the identity mapping, shows how the coordinate system of $\mu$ should be deformed to be registered to the coordinate system of $\mu^{i}$.

This raises the question of how to characterise the optimal transportation maps. For instance, in the one-dimensional case, if $\mu$ and $\nu$ are probability measures on $\mathbb{R}$, and $\mu$ is diffuse we may write

[TABLE]

where $G_{\mu}(t)={\int_{-\infty}^{t}\!\,\mathrm{d}\mu(x)}$, $G_{\nu}(t)={\int_{-\infty}^{t}\!\,\mathrm{d}\nu(x)}$ are their distribution functions and $G^{-1}_{\nu}$ is the quantile function of $\nu$. This characterises optimal maps in one dimension as non-decreasing functions. More generally, when one has measures on $\mathbb{R}^{d}$, the class of optimal maps can be seen to be that of monotone maps (see Section 5.5), defined as fields $\mathbf{t}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ that are obtained as gradients of convex functions $\varphi:\mathbb{R}^{d}\rightarrow\mathbb{R}$,

[TABLE]

This is known as Brenier’s characterisation (Villani [74, Theorem 2.12]). With these basic definitions in place, we are now ready to consider the problem of finding a Fréchet mean of a collection of measures – the latter viewed as the common template measure that was deformed to give rise to these measures.

3 Sample Setting

3.1 Fréchet Means and Optimal Registration

The notion of a Fréchet mean (Fréchet [32]) generalises that of the mean in a normed vector space to a general metric space. Though it has primarily been studied on Riemannian manifolds, the generality of its definition allows it to be used very broadly: it replaces the usual “sum of squares”, with a “sum of squared distances”, the Fréchet functional. A closely related notion is that of a Karcher mean (Karcher [45]; Le [50]), a term that describes stationary points of the sum of squares functional, when the latter is differentiable. See Kendall [46], and Kendall & Le [47] for an overview and a detailed review, respectively. In the context of Wasserstein space, a Fréchet mean of a collection of measures $\{\mu^{1},\ldots,\mu^{N}\}$, is a minimiser of the Fréchet functional

[TABLE]

over elements $\gamma$ in the Wasserstein space $\mathcal{P}_{2}(\mathbb{R}^{d})$, and a Karcher mean is a stationary point of $F$. The functional will be finite for any $\gamma\in\mathcal{P}_{2}(\mathbb{R}^{d})$, provided that it is so for some $\gamma_{0}$. Population versions, assuming $\mathcal{P}_{2}(\mathbb{R}^{d})$ is endowed with a probability measure, can also be defined, replacing summation by expectation with respect to that law. Interestingly, Fréchet himself [33] considered the Wasserstein metric between probability measures on $\mathbb{R}$, and some refer to this as the Fréchet distance (e.g. Dowson & Landau [29]). In general, existence and uniqueness of a sample Fréchet mean can be subtle, but Agueh & Carlier [2] have shown that it will uniquely exist in the Wasserstein space, provided that some regularity is asserted222For a population version, one needs to tackle measurability and identifiability issues, see Section 4.1. Here and in the following, we call a measure regular if it is absolutely continuous with respect to Lebesgue measure (this condition can be slightly weakened [2]).

Proposition 1** (Agueh & Carlier [2]).**

Let $\{\mu^{1},\ldots,\mu^{N}\}$ be a collection in the Wasserstein space of measures $\mathcal{P}_{2}(\mathbb{R}^{d})$. If at least one of the measures is regular with bounded density, then their Fréchet mean exists, is unique, and is regular.

We will now show that, once the Fréchet mean $\bar{\mu}$ of $\{\mu^{1},\ldots,\mu^{N}\}$ has been determined, it may be used to optimally multi-couple the measures $\{\mu^{1},\ldots,\mu^{n}\}$ in $\mathbb{R}^{d\times N}$, in terms of pairwise mean square distances, thus providing a solution to the multidimensional Monge–Kantorovich problem considered by Gangbo & Świȩch [36]. That is, $\bar{\mu}$ can be used to construct a random vector whose marginals are as concentrated as possible in terms of pairwise mean-square distance, subject to the constraint of having laws $\{\mu^{1},\ldots,\mu^{N}\}$.

Our first result combines results of [2] and [36] to illustrate precisely how (also see Pass [62, Theorem 4.2.2] for an analogous result when considering continuous flows of measures).

Proposition 2** (Optimal Multicoupling via Fréchet Means).**

Let $\{\mu^{1},\ldots,\mu^{N}\}$ be regular probability measures in $\mathcal{P}_{2}(\mathbb{R}^{d})$, one with bounded density, and let $\bar{\mu}$ be their (unique) Fréchet mean with respect to the Wasserstein metric. Let $Z\sim\bar{\mu}$ and define

[TABLE]

where $\mathbf{t}_{\bar{\mu}}^{\mu^{i}}$ is the optimal transport plan pushing $\bar{\mu}$ forward to $\mu^{i}$. Then $X_{i}\sim\mu^{i}$ for $i=1,\ldots,N$ and furthermore,

[TABLE]

for any other $\bm{Y}=(Y_{1},\ldots,Y_{N})$ such that $Y_{i}\sim\mu^{i}$, $i=1,\ldots,N$.

In the language of shape theory, the Fréchet mean $\bar{\mu}$ may be used as a template to jointly register the collection of measures, just as Euclidean configurations can be registered to their Procrustes mean by a Procrustes analysis (Goodall [37]). Only in this case, instead of the similarity group of shape theory, registration is deformation based, by means of the collection of maps $\{\mathbf{t}_{\bar{\mu}}^{\mu^{i}}\}_{i=1}^{N}$, where $\mathbf{t}_{\bar{\mu}}^{\mu^{i}}$ is the optimal transport map

[TABLE]

By analogy to shape theory, we shall refer to these as Procrustes maps. These yield a common coordinate system (corresponding to $\bar{\mu}$) where one can best compare samples from each measure, similarly to “quantile renormalisation” in one dimension, e.g. Bolstad et al. [17], Gallon et al. [35]. The Procrustes maps can also be used in order to produce a Principal Component Analysis, capturing the main modes of deformation-based variation (Bigot et al. [13], Panaretos & Zemel [61]; Huckemann, Munk & Hotz [42], Wang et al. [75]).

3.2 Wasserstein Geometry and the Gradient of the Fréchet Functional

In this section, we determine the conditions for the Fréchet derivative of the Fréchet functional (3.1) to be well defined, and determine its functional form. Furthermore, we characterise Karcher means and give criteria for their optimality, opening the way for the determination of the Fréchet mean. The key to our analysis will be to exploit the tangent bundle over the Wasserstein space of regular measures.

3.2.1 The Tangent Bundle

Let $\mathcal{P}_{2}(\mathbb{R}^{d})$ be the Wasserstein space of probability measures $\mu$ on $\mathbb{R}^{d}$ such that ${\int_{\mathbb{R}^{d}}\!\|x\|^{2}\,\mathrm{d}\mu(x)}$ is finite, as defined in Section 2. An absolutely continuous measure on $\mathbb{R}^{d}$ will be called regular. When $\mu^{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ is regular and $\mu^{1}\in\mathcal{P}_{2}(\mathbb{R}^{d})$, the transportation map $\mathbf{t}_{\mu^{0}}^{\mu^{1}}$ uniquely exists, in which case there is a unique geodesic curve between $\mu^{0}$ and $\mu^{1}$. Using again the notation $\mathbf{i}$ for the identity map, this geodesic is given by

[TABLE]

This curve is known as McCann’s interpolation (McCann [54], Villani [74]). The tangent space at an arbitrary $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ is then (Ambrosio et al. [7, Definition 8.4.1, p. 189])

[TABLE]

where $C_{c}^{\infty}(\mathbb{R}^{d})$ denotes infinitely differentiable functions $\varphi:\mathbb{R}^{d}\to\mathbb{R}$ with compact support, and the closure operation is taken with respect to the space $L^{2}(\mu)$. Note the interesting fact that the closure operation is the only aspect of the tangent space that directly involves the measure $\mu$. An equivalent definition, which is more useful to us, is given by Ambrosio et al. [7, Definition 8.5.1, p. 195]:

[TABLE]

that is, we take the collection of $\mathbf{r}$’s that are optimal maps from $\mu$ to $\mathbf{r}\#\mu$; i.e. the gradients of convex functions. This is a linear space (not just a cone) by the first definition, even though it is not obvious from the second. The definitions are equivalent by Theorem 8.5.1 of Ambrosio et al. [7, p. 195]. As was mentioned above, when $\mu^{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ is regular, every measure $\mu^{1}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ admits a unique optimal map $\mathbf{t}_{\mu^{0}}^{\mu^{1}}$ that pushes $\mu^{0}$ forward to $\mu^{1}$. Thus, the exponential map

[TABLE]

is surjective, and its inverse, the log map

[TABLE]

is well-defined throughout $\mathcal{P}_{2}(\mathbb{R}^{d})$. In particular, the geodesic $\left[\mathbf{i}+t(\mathbf{t}_{\mu^{0}}^{\mu^{1}}-\mathbf{i})\right]\#\mu^{0}$ is mapped bijectively to the line segment $t(\mathbf{t}_{\mu^{0}}^{\mu^{1}}-\mathbf{i})\in\mathrm{Tan}_{\mu^{0}}$ through the log map.

3.2.2 Gradient of the Fréchet functional

We will now exploit the tangent bundle structure described in the previous section in order to determine the gradient of the empirical Fréchet functional. Fix $\mu^{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ and consider the function

[TABLE]

When $\mu$ is regular, we have that ([7, Corollary 10.2.7, p. 239]), for any $\mu^{0}$

[TABLE]

where the convergence $\nu\to\mu$ is with respect to the Wasserstein distance. The integral above can be seen as the inner product

[TABLE]

in the space $L^{2}(\mu)$ that includes as a (closed) subspace the tangent space $\mathrm{Tan}_{\mu}$. In terms of this inner product and the log map, we can write

[TABLE]

so that $F_{0}$ is Fréchet-differentiable at $\mu$ with derivative

[TABLE]

We have proven:

Theorem 1** (Gradient of the Fréchet Functional).**

Fix a collection of measures $\mu^{1},\dots,\mu^{N}\in\mathcal{P}_{2}(\mathbb{R}^{d})$. When $\gamma$ is regular, the Fréchet functional

[TABLE]

is Fréchet-differentiable, and its gradient satisfies

[TABLE]

3.2.3 Karcher and Fréchet Means

We can now characterise Karcher means, and also show that the empirical Fréchet mean must be sought amongst them, by an immediate corollary to Theorem 1:

Corollary 1**.**

Let $\mu^{1},\dots,\mu^{N}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be regular measures, one of which with bounded density. A measure $\mu$ is a Karcher mean of $\{\mu^{i}\}$ if and only if

[TABLE]

Furthermore, the Fréchet mean of $\{\mu^{i}\}$ is itself a Karcher mean.

In fact, the corollary suggests that a Karcher mean is “almost” a Fréchet mean: Agueh and Carlier [2] show by convex optimisation methods that if $\sum_{i=1}^{N}\left(\mathbf{t}_{\mu}^{\mu^{i}}-\mathbf{i}\right)=0$ everywhere on $\mathbb{R}^{d}$ (rather than just $\mu$-almost everywhere), then $\mu$ is in fact the unique Fréchet mean. Thus one hopes that this “gap of measure zero” can be bridged: that a sufficiently regular Karcher mean should in fact be a Fréchet mean. We now show that this is indeed the case; if $\mu^{1},\dots,\mu^{N}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ are smooth measures with convex support, then a smooth Karcher mean of same support must be the unique Fréchet mean:

Theorem 2** (Optimality Criterion for Karcher Means).**

Let $\mu^{i}$ for $i=1,\dots,N$ be probability measures on an open convex $X\subseteq\mathbb{R}^{d}$ whose densities $g^{i}$ are bounded and strictly positive on $X$ and let $\mu$ be a regular Karcher mean of $\{\mu^{i}\}$ with density $f$. Then $\mu$ is the unique Fréchet mean of $\{\mu^{i}\}$, provided one of the following holds:

$X=\mathbb{R}^{d}$, $f$ is bounded and strictly positive, and the densities $f,g^{1},\dots,g^{N}$ are of class $C^{1}$; 2. 2.

$X$* is bounded, $\mu(X)=1$, $f$ is bounded, and the densities $f,g^{1},\dots,g^{N}$ are bounded from below on $X$.*

Remark 1**.**

In the first condition, the $C^{1}$ assumption can be weakened to Hölder continuity of the densities for some exponent $\alpha\in(0,1]$.

Remark 2**.**

We conjecture that a stronger result should be valid: specifically, if $\mu^{1},\dots,\mu^{N}$ satisfy the conditions of Theorem 2, then we conjecture the Fréchet functional $F$ to in fact have a unique Karcher mean, coinciding with the Fréchet mean.

3.3 Gradient Descent and Procrustes Analysis

3.3.1 Elements of the Algorithm

Let $\mu^{1},\dots,\mu^{N}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be regular and let $\gamma_{j}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be a regular measure, representing our current estimate of the Fréchet mean of $\mu^{1},\dots,\mu^{N}$ at step $j$. Following the discussion above, it makes sense to introduce a step size $\tau_{j}>0$, and to carry out a steepest descent in the space of measures (e.g. Molchanov & Zuyev [55]), following the negative of the gradient:

[TABLE]

In order to guarantee that the descent is well-defined, we must make sure that the gradient itself will remain well-defined as we iterate over $j$. In view of Theorem 1, this requires showing that $\gamma_{j+1}$ remains regular whenever $\gamma_{j}$ is regular. This is indeed the case, at least if the step size is contained in $[0,1]$:

Lemma 1** (Regularity of the iterates).**

If $\gamma_{0}$ is regular and $\tau_{0}\in[0,1]$ then so is $\gamma_{1}$.

Lemma 1 suggests that the step size must be restricted to $[0,1]$. The next result suggests that the objective function essentially tells us that the optimal step size, achieving the maximal reduction of the objective function (thus corresponding to an approximate line search), is exactly equal to 1:

Lemma 2** (Optimal Stepsize).**

If $\gamma_{0}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ is regular then

[TABLE]

and the bound on the right-hand side of the last display is minimised when $\tau=1$.

In light of the results in Lemmas 1 and 2, one needs only concentrate on the case $\tau_{j}=1$. This has an interesting ramification: when $\tau=1$, the gradient descent iteration is structurally equivalent to a Procrustes analysis. Specifically, the gradient descent algorithm proceeds by iterating the two steps of a Procrustes analysis (Gower [38]; Dryden & Mardia [30, p. 90]):

1. (1)

Registration: Each of the measures $\{\mu^{1},\ldots,\mu^{N}\}$ is registered to the current template $\gamma_{j}$, via the optimal transportation (registration) maps $\mathbf{t}_{\gamma_{j}}^{\mu^{i}}$. In geometrical terms, the measures $\{\mu^{1},\ldots,\mu^{N}\}$ are lifted to the tangent space at $\gamma_{j}$ (via the log map), and their linear representation on the tangent space is expressed in local coordinates which coincide with the maps $\mathbf{t}_{\gamma_{j}}^{\mu^{i}}-\mathbf{i}=\log_{\gamma_{j}}(\mu^{i})$. These can be seen as a common coordinate system for $\{\mu^{1},\ldots,\mu^{N}\}$, i.e. a registration. 2. (2)

Averaging: The registered measures are averaged coordinate-wise, using the common coordinates system by the registration step (1). In geometrical terms, the linear representation of $\{\mu^{1},\ldots,\mu^{N}\}$ afforded by their local coordinates $\mathbf{t}_{\gamma_{j}}^{\mu^{i}}-\mathbf{i}=\log_{\gamma_{j}}(\mu^{i})$ is averaged linearly. The linear average is then retracted back onto the manifold via the exponential map to yield the estimate at the $(j+1)$-step.

That the gradient descent reduces to Procrustes analysis is not simply of aesthetic value. It is of the essence, as it shows that the algorithm relies entirely on solving a succession of pairwise optimal transportation problems, thus reducing the determination of the Fréchet mean to the classical Monge problem of optimal transportation (e.g. Benamou and Brenier [10], Haber et al. [40], Chartrand et al. [23]). After all, this is precisely the point of a Procrustes algorithm: exploiting the (easier) problem of pairwise registration to solve the (harder problem) of multi-registration. We note that, further to requiring the ability to solve the pairwise optimal transportation problem, and the regularity conditions on the measures, the algorithm does not require additional structural assumptions/workarounds to reduce the problem to the one-dimensional case (as in, for example the “admissibility” approach of Boissard et al. [16]). An additional practical advantage is that Procrustes algorithms are easily parallelisable, since one can distribute the solution of the pairwise transport problems at each step $j$. Any regular measure can serve as an initial point for the algorithm, for instance one of the $\mu^{i}$. We should mention at this point that, if one is content with obtaining an approximate or regularised Fréchet mean, then there are several numerical strategies available, and there is a rapidly growing literature for the efficient computation of such schemes – we briefly summarise some such approaches in the concluding remarks section (Section 7).

The gradient/Procrustes iteration is presented succinctly as Algorithm 1.

3.3.2 Convergence of the Algorithm

In order to tackle the issue of convergence, we will use an approach that is specific to the nature of optimal transportation. The reason is that Hessian type arguments that are used to prove similar convergence results for gradient descent on Riemmanian manifolds (Afsari et al. [1]) or Procrustes algorithms (Le [52], Groisser [39]) do not apply here, since the Fréchet functional may very well fail to be twice differentiable. Still, this specific geometry of Wasserstein space affords some advantages; for instance, we will place no restriction on the starting point for the iteration, except that it be regular:

Theorem 3** (Limit Points are Karcher Means).**

Let $\mu^{1},\dots,\mu^{N}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be absolutely continuous probability measures, one of which with bounded density. Then, the sequence generated by Algorithm 1 stays in a compact set of the Wasserstein space $\mathcal{P}_{2}(\mathbb{R}^{d})$, and any limit point of the sequence is a Karcher mean of $(\mu^{1},\dots,\mu^{N})$.

In view of Corollary 1, this immediately implies:

Corollary 2** (Wasserstein Convergence of Gradient Descent).**

Under the conditions of Theorem 3, if $F$ has a unique stationary point, then the sequence $\{\gamma_{j}\}$ generated by Algorithm 1 converges to the Fréchet mean of $\{\mu^{1},\dots,\mu^{N}\}$ in the Wasserstein metric,

[TABLE]

Of course, combining Theorem 3 with Theorem 2 shows that the conclusion of Corollary 2 holds when the appropriate assumptions on $\{\mu^{i}\}$ and the Karcher mean $\mu$ are satisfied. The proof of Theorem 3 is elaborate, and is constructed via a series of intermediate results in a separate section (Section 5.3.1) in the interest of tidiness. The main challenge is that the standard condition used for convergence of gradient descent algorithms, that gradients be Lipschitz, fails to hold in this setup. Indeed, $F$ is not differentiable on discrete measures, and these constitute a dense subset of the Wasserstein space.

3.3.3 Uniform Convergence of Procrustes Maps and Multicoupling

We conclude our analysis of the algorithm by turning to the Procrustes maps $\mathbf{t}^{\bar{\mu}}_{\mu^{i}}$, which optimally couple each sample observation $\mu^{i}$ to their Fréchet mean $\bar{\mu}$. These are the key objects required for the solution of the multicoupling problem (as established in Proposition 2), and one would use the limit of $\mathbf{t}^{\gamma_{j}}_{\mu^{i}}$ in $j$ as their approximation. However, the fact that $d(\gamma_{j},\bar{\mu}){\rightarrow}0$ does not immediately imply the convergence of $\mathbf{t}^{\gamma_{j}}_{\mu^{i}}$ to $\mathbf{t}^{\bar{\mu}}_{\mu^{i}}$: the Wasserstein convergence only means that certain integrals of the warp maps converge. Still, convergence of the warp maps does hold, indeed uniformly so on compacta, $\bar{\mu}$-almost everywhere:

Theorem 4** (Uniform Convergence of Procrustes Maps).**

Under the conditions of Corollary 2, there exist sets $A,B^{1},\dots,B^{N}\subseteq\mathbb{R}^{d}$ such that $\bar{\mu}(A)=1=\mu^{1}(B^{1})=\dots=\mu^{N}(B^{N})$ and

[TABLE]

for any pair of compacta $\Omega_{1}\subseteq A$, $\Omega_{2}\subseteq B^{i}$, where the sequence $\mathbf{t}^{\gamma_{j}}_{\mu^{i}}$ and $\mathbf{t}_{\gamma_{j}}^{\mu^{i}}=\left(\mathbf{t}^{\gamma_{j}}_{\mu^{i}}\right)^{-1}$ are the Procrustes maps generated by Algorithm 1. If in addition all the measures $\{\mu^{1},\dots,\mu^{N}\}$ have the same support, then one can choose the sets so that $B^{1}=\dots=B^{N}$.

With both ingredients of the registration problem in hand, we deduce a solution to the latter:

Corollary 3** (Convergence of Multicouplings).**

Under the conditions of Corollary 2, the sequence of multicouplings

[TABLE]

of $\{\mu^{1},\dots,\mu^{N}\}$ converges (in Wasserstein distance on $(\mathbb{R}^{d})^{N}$) to the optimal multicoupling $(\mathbf{t}_{\bar{\mu}}^{\mu^{1}},\dots\mathbf{t}_{\bar{\mu}}^{\mu^{n}})\#\bar{\mu}$.

4 Population Setting

In order to carry out inference, we must relate the sample collection of measures to a population, and show that the relevant quantities are identifiable parameters. Furthermore, in practice the sample measures will only be discretely observed, and this must be taken into account. We now formulate such a model, and study its nonparametric estimation from discrete observations.

4.1 Deformation Models and Discrete Observation

Let $\lambda$ be a regular probability measure with a strictly positive density on a convex compact $K\subset\mathbb{R}^{d}$ of positive Lebesgue measure333In applied settings, the point processes will be observed on a bounded observation window $K$. For this reason as well as the sake of simplicity, we restrict our discussion to a given compact set (but remark that it could be extended to unbounded observation windows subject to further conditions). , and let $\{\Pi_{1},\ldots,\Pi_{N}\}$ be i.i.d point processes with intensity measure $\lambda$,

[TABLE]

for all Borel subsets $A\subseteq K$. Instead of observing the true processes $\{\Pi_{1},\ldots,\Pi_{N}\}$, we are able to observe warped versions

[TABLE]

with conditional warped mean measures

[TABLE]

where the $\{T_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\}$ are i.i.d random homeomorphisms on $K$, satisfying the properties of

Unbiasedness: the Fréchet mean of $\Lambda_{i}=T_{i}\#\lambda$ is $\lambda$. 2. 2.

Regularity: $T_{i}$ is a gradient of a convex function on $K$.

The conditional mean measures $\{\Lambda_{i}=T_{i}\#\lambda\}_{i=1}^{N}$ play the role of the unobservable sample of random measures generated from a population law constructed via random deformations of the template $\lambda$. The processes $\{\widetilde{\Pi}_{i}\}_{i=1}^{N}$ play the role of the discretely observed versions of the $\{\Lambda_{i}\}_{i=1}^{N}$. Conditions (1) and (2) state that the deformations $\{T_{i}\}$ are identifiable. They can also be motivated from first principles: (1) states that the maps do not deform the template $\lambda$ on average (otherwise this “average deformation” would be by definition the template); and (2) states that among all possible deformations that could have mapped $\lambda$ to $\Lambda_{i}$, we take the parsimonious choice of the optimal deformation. The importance and canonicity of these two assumptions has been discussed in depth in Panaretos & Zemel [61, Section 3.3], who study a one-dimensional version of the above problems (which is qualitatively very different, given the flat nature of 1d Wasserstein space, and the availability of explicit closed form expressions).

The connection of this deformation model to Fréchet means, via the optimal maps, is now given as follows (in a general setup, encompassing our model setup). Let $C_{b}(K,\mathbb{R}^{d})$ be the space of continuous bounded functions $f:K\to\mathbb{R}^{d}$ endowed with the supremum norm $\|f\|_{\infty}=\sup_{x\in K}\|f(x)\|$.

Theorem 5** (mean identity warp functions and Fréchet means).**

Let $K\subset\mathbb{R}^{d}$ be a compact convex set of positive Lebesgue measure, and let $\lambda\in\mathcal{P}_{2}(K)$ be regular. Consider the random measure $\Lambda=T\#\lambda$, where $T:K\to K$ is a random deformation (viewed as a random element in $C_{b}(K,\mathbb{R}^{d})$), almost surely injective, and satisfying

almost surely there exists a convex function $\phi$ such that $T=\nabla\phi$ on the interior of $K$; 2. 2.

$\mathbb{E}[T(x)]=x$* for all $x\in K$ (or on a dense subset of $K$);* 3. 3.

almost surely $T$ is differentiable with a nonsingular derivative for almost all $x\in K$.

Then $\lambda$ is the unique Fréchet mean of $\Lambda$, i.e., the unique minimiser of the population Fréchet functional $\gamma\mapsto\mathbb{E}d^{2}(\Lambda,\gamma)$.

An important requirement for the statement and proof of Theorem 5 is that $\phi$, $\phi^{*}$ and $\Lambda$ are measurable as random elements in the appropriate spaces; this is not a priori obvious, but is established as part of the proof.

The statistical problem will now be to estimate the unknown structural mean measure $\lambda$, and the registration maps $T_{i}$ non-parametrically, by smoothing the observed point processes $\{\widetilde{\Pi}_{1},\ldots,\widetilde{\Pi}_{N}\}$. Once $\lambda$ and $\{T_{i}\}$ have been estimated, the processes $\{\widetilde{\Pi}_{1},\ldots,\widetilde{\Pi}_{N}\}$ can be registered by applying the inverses of the estimated maps $T_{i}$, allowing for further analysis of the point processes in a functional data context. Theorem 5 guarantees that the estimands considered are identifiable.

4.2 Regularised Nonparametric Estimation

In order to estimate the $\lambda$ and the $\{\Lambda_{i},T_{i}\}$, we will follow the steps below:

Regularisation: Estimate $\Lambda_{i}=T_{i}\#\lambda$ by a regular kernel estimator $\widehat{\Lambda}_{i}$ restricted on $K$,

[TABLE]

where $\psi:\mathbb{R}^{d}\to(0,\infty)$ is a unit-variance isotropic density function, $\psi_{\sigma}(x)=\sigma^{-d}\psi(x/\sigma)$ for $\sigma>0$ (more generally, $\psi$ could be non-isotropic, having a bandwidth matrix, but we focus on the isotropic case for simplicity), and $\widetilde{\Pi}_{i}$ is the sum of dirac masses $\sum_{j=1}^{m}\delta\{x_{i}\}$. If $\widetilde{\Pi}_{i}$ contains no points (that is, $m=0$), define $\hat{\Lambda}_{i}$ to be the (normalised) Lebesgue measure on $K$. 2. 2.

Fréchet Mean Estimation: Estimate $\lambda$ by the empirical Fréchet mean $\hat{\lambda}$ of $\widehat{\Lambda}_{1},\dots,\widehat{\Lambda}_{N}$, using the Procrustes Algorithm 1. 3. 3.

Procrustes Analysis: Estimate $T_{i}$ by the optimal transportation map of $\widehat{\lambda}$ onto $\widehat{\Lambda}_{i}$, as given by the final step in the iteration of Algorithm 1. Estimate the map $T_{i}^{-1}$ by $\widehat{T^{-1}_{i}}=\widehat{T}^{-1}_{i}$. 4. 4.

Registration: Register the observed point processes to a common coordinate system by defining $\widehat{\Pi}_{i}=\widehat{{T^{-1}_{i}}}{\#}\widetilde{\Pi}_{i}$.

In the next section, we will prove that our estimates are consistent for their population version, as the number of observed processes, and the number of points per process diverge.

4.3 Asymptotic Theory

To establish consistency, we will use the dense asymptotics regime of functional data analysis, adapted to the current setting. We will consider a setup where the number of observed point processes $n$ diverges, and the (mean) number of points in each observed process, $\mathbb{E}[\widetilde{\Pi}_{i}(K)]$, diverge too. Here we use the index notation “$n$” rather than “$N$” to emphasize that the index is no longer held fixed. Specifically, let $(\Pi_{1}^{(n)},\Pi_{2}^{(n)},\dots,\Pi_{n}^{(n)})_{n=1}^{\infty}$ be a triangular array of row-independent and identically distributed point processes on $K$ following the same infinitely divisible distribution and having mean measure $\tau_{n}\lambda$, where $\tau_{n}>0$ are constants. Let $T_{1},\dots,T_{n}$ be independent and identically distributed realisations of a random homeomorphism $T$ of $K$ satisfying the unbiasedness and regularity assumptions of Section 4.1. Let $\widetilde{\Pi}_{i}^{(n)}=T_{i}\#\Pi_{i}^{(n)}$ and set $\Lambda_{i}=T_{i}\#\lambda=\tau_{n}^{-1}\mathbb{E}[\widetilde{\Pi}_{i}^{(n)}|T_{i}]$. Suppose that $\widehat{\Lambda}_{i}$ is an estimator of $\Lambda_{i}$, constructed by kernel smoothing of $\Pi_{i}^{(n)}$ using a (possibly random) bandwidth $\sigma_{i}^{(n)}$, as described in the previous section. Correspondingly, let $\widetilde{\Pi}_{i}^{(N)}=T_{i}\#\Pi_{i}^{(n)}$ and set $\Lambda_{i}=T_{i}\#\lambda=\tau_{n}^{-1}\mathbb{E}[\widetilde{\Pi}_{i}^{(n)}|T_{i}]$.

Theorem 6** (Consistency of the regularised Fréchet Mean).**

If $\tau_{n}/\log n\to\infty$ and $\sigma_{n}=\max_{i}\sigma^{(n)}_{i}\stackrel{{\scriptstyle\rm{p}}}{{\to}}0$ then

For any $i$,

[TABLE] 2. 2.

The estimator $\widehat{\lambda}_{n}$ is strongly consistent

[TABLE]

If the smoothing is carried out independently across trains, that is, $\sigma_{i}^{(n)}$ depends only on $\widetilde{\Pi}_{i}^{(n)}$, then the result still holds if merely $\tau_{n}\to\infty$.

If $\mathbb{E}\left[\Pi_{1}^{(1)}\right]^{4}<\infty$, $\sum_{n}\tau_{n}^{-2}<\infty$ and $\sigma_{n}\stackrel{{\scriptstyle\rm{as}}}{{\to}}0$ then convergence almost surely holds.

Remark 3**.**

There is no lower bound on $\sigma_{n}$, and it can vanish at any rate, provided it is strictly positive. In practice, however, if $\sigma_{n}$ is very small, then the densities of $\widehat{\Lambda}_{i}$ will have very high peaks, and the constant $C_{\mu}$ in Proposition 4 (with $\mu^{i}=\widehat{\Lambda}_{i}$) will be large (essentially proportional to $1/\sigma_{n}$). The proof of Proposition 3 suggests that this may slow down the convergence of Algorithm 1.

Remark 4**.**

It is worth remarking that Le Gouic & Loubes [53, Theorem 3] consider the stability of Fréchet means in a rather general setting; verification of their assumptions in our particular setting, however, is quite involved and in fact essentially amounts to directly proving Theorem 6.

Our next two results concern the (uniform) consistency of the Procrustes registration procedure. Though the results themselves parallel their one-dimensional counterparts (see Panaretos & Zemel [61]), their proofs are entirely different, and substantially more involved (because the geometry of monotone mappings in $\mathbb{R}^{d}$ is far more rich than the geometry of monotone maps on $\mathbb{R}$). In particular, we have:

Theorem 7** (Consistency of Procrustes Maps).**

Under the same conditions of Theorem 6, for any $i$ and any compact set $\Omega\subseteq\mathrm{int}(K)$,

[TABLE]

The same remarks at the end of the statement of Theorem 6 apply here as well.

Corollary 4** (Consistency of Procrustes Registration).**

Under the same conditions of Theorem 6, the registration procedure is consistent: for any $i$

[TABLE]

provided one of the following conditions holds:

Every point of the boundary of $K$ is exposed, that is, for any $y\in\partial K$ there exists $\alpha\in\mathbb{R}^{d}$ such that

[TABLE] 2. 2.

The warp map $T_{i}$ is strictly monotone

[TABLE]

The first condition is satisfied by any ellipsoid in $\mathbb{R}^{d}$ and more generally if the boundary of $K$ can be written as $\partial K=\{x:\varphi_{K}(x)=0\},$ for a strictly convex function $\varphi_{K}$. Indeed, if $\alpha$ creates a supporting hyperplane to $K$ at $y$ and $\langle\alpha,y\rangle=\langle\alpha,y^{\prime}\rangle$ for $y\neq y^{\prime}$, then as $\varphi_{K}$ is strictly convex on the line segment $[y,y^{\prime}]$, it is impossible that $y^{\prime}\in K$ without the hyperplane intersecting the interior of $K$. Although this condition excludes some interesting cases, perhaps most prominently polyhedral sets such as $K=[0,1]^{d}$, such sets can be approximated by convex sets that do satisfy it (Krantz [49, Proposition 1.12]).

As for the second condition, in general it will hold almost surely. Indeed, as $T_{i}\#\lambda=\Lambda_{i}$ and both measures are absolutely continuous, there exists a $\lambda$-null set $\mathcal{N}$ such that $T_{i}$ is strictly monotone outside $\mathcal{N}$ [7, Proposition 6.2.12]. By assumption $\lambda$ has a strictly positive density on $K$, so that $\lambda$-null subsets of $K$ are precisely the Lebesgue null subsets of $K$. In that sense, this condition is not overly restrictive, and will most likely be satisfied under additional regularity assumptions on the warp maps $T_{i}$ and, possibly, $K$.

5 Proofs of Formal Statements

Our proofs will require us to establish some analytical results that are intrinsic to the optimal transportation problem. These are essential for the proofs, especially of our main results, and some are non-trivial. For tidiness, we will state and prove these results separately at the end of this section (Section 5.5), developing our main results first, and referring to the analytical background when necessary.

5.1 Proofs of Statements in Section 3.1

Proof of Proposition 2.

The optimisation problem

[TABLE]

is equivalent to minimising

[TABLE]

and Agueh & Carlier [2, Proposition 4.2] show that $\min_{\mu}F(\mu)=\min_{\xi}G(\xi)$.

Since $\bar{\mu}$ is regular [2, Proposition 5.1], $\bm{X}$ is well-defined and has joint distribution

[TABLE]

Since the coordinates of $h$ have mean identity (see [2, Equation (3.9)] or Corollary 1),

[TABLE]

Thus $\xi^{\prime}$ is optimal. ∎

5.2 Proofs of Statements in Section 3.2

Proof of Corollary 1.

The characterisation of Karcher means is immediate from Theorem 1. The fact that the Fréchet mean $\mu$ satisfies $\sum_{i=1}^{N}\left(\mathbf{t}_{\mu}^{\mu^{i}}-\mathbf{i}\right)=0$ $\mu$-almost everywhere follows by a result of Agueh & Carlier [2]. For an alternative proof using the tangent bundle, see the supplementary material (Section 8). ∎

Proof of Theorem 2.

The result exploits Caffarelli’s regularity theory for Monge–Ampère equations. In the first case, by Theorem 4.14(iii) in Villani [74] there exist $C^{1}$ (in fact, $C^{2,\alpha}$) convex potentials $\varphi_{i}$ on $\mathbb{R}^{d}$ with $\mathbf{t}_{\mu}^{\mu^{i}}=\nabla\varphi_{i}$, so that $\mathbf{t}_{\mu}^{\mu^{i}}(x)$ is a singleton for all $x\in\mathbb{R}^{d}$. The set $\{x\in\mathbb{R}^{d}:\sum\mathbf{t}_{\mu}^{\mu^{i}}(x)/N\neq x\}$ is $\mu$-negligible (and hence Lebesgue-negligible) and open by continuity. It is therefore empty, so $F^{\prime}(\mu)=0$ everywhere, and $\mu$ is the Fréchet mean (see the discussion after Corollary 1).

In the second case, by the main theorem in Caffarelli [21, p. 99], and the same argument, we have $\sum\mathbf{t}_{\mu}^{\mu^{i}}(x)/N=x$ for all $x\in X$. Since $X$ is convex, there must exist a constant $C$ such that $\sum\varphi_{i}(x)=C+N\|x\|^{2}/2$ for all $x\in X$. Hence Equation (3.9) in [2] holds with $\mathbb{R}^{d}$ replaced by $X$. Repeating the proof of Proposition 3.8 in [2], we see that $\mu$ minimises $F$ on $\mathcal{P}_{2}(X)$, the set of measures supported on $X$. (All the integrals that appear in the proof can be taken on $X$, where we know the inequality holds). Again by convexity of $X$, the minimiser of $F$ must be444We know that the minimiser must be in $\mathcal{P}_{2}(\overline{X})$, but minimising on $\mathcal{P}_{2}(X)$ suffices by continuity of $F$. in $\mathcal{P}_{2}(X)$ (see the existence proof at the beginning of the proof of Theorem 5 in the supplementary material, Section 8). ∎

5.3 Proofs of Statements in Section 3.3

Proof of Lemma 1.

By [7, Proposition 6.2.12] there exists a $\gamma_{0}$-null set $A_{i}$ such that on $\mathbb{R}^{d}\setminus A_{i}$, $\mathbf{t}_{\gamma_{0}}^{\mu^{i}}$ is differentiable, $\nabla\mathbf{t}_{\gamma_{0}}^{\mu^{i}}>0$ (positive definite), and $\mathbf{t}_{\gamma_{0}}^{\mu^{i}}$ is strictly monotone

[TABLE]

Since $\mathbf{t}_{\gamma_{0}}^{\gamma_{1}}=(1-\tau)\mathbf{i}+\tau N^{-1}\sum_{i=1}^{N}\mathbf{t}_{\gamma_{0}}^{\mu^{i}}$, it stays strictly monotone (hence injective) and $\nabla\mathbf{t}_{\gamma_{0}}^{\gamma_{1}}>0$ outside $A=\cup A_{i}$, which is a $\gamma_{0}$-null set.

Let $h_{0}$ denote the density of $\gamma_{0}$ and set $\Sigma=\mathbb{R}^{d}\setminus A$. Then $\mathbf{t}_{\gamma_{0}}^{\gamma_{1}}|_{\Sigma}$ is injective and $\{h_{0}>0\}\setminus\Sigma$ is Lebesgue negligible because

[TABLE]

and the integrand is strictly positive. Since $|\mathrm{det}\nabla\mathbf{t}_{\gamma_{0}}^{\mu^{i}}|>0$ on $\Sigma$ we obtain that $\gamma_{1}=\mathbf{t}_{\gamma_{0}}^{\mu^{i}}\#\gamma_{0}$ is absolutely continuous by [7, Lemma 5.5.3]. ∎

Proof of Lemma 2.

Let $S_{i}=\mathbf{t}_{\gamma_{0}}^{\mu^{i}}$ be the optimal map from $\gamma_{0}$ to $\mu^{i}$, and set $W_{i}=S_{i}-\mathbf{i}$. Then

[TABLE]

with the inner product being in $L^{2}(\gamma_{0})$. By definition

[TABLE]

This is a map that pushes forward $\mu^{i}$ to $\gamma_{1}$ (not necessarily optimally). Hence

[TABLE]

Now $\mu^{i}=S_{i}\#\gamma_{0}$, which means that ${\int\!f\,\mathrm{d}\mu^{i}}={\int\!(f\circ S_{i})\,\mathrm{d}\gamma_{0}}$ for any measurable $f$. This change of variables gives

[TABLE]

The norm is always in $L^{2}(\gamma_{0})$, regardless of $i$. Developing the squares, summing over $i=1,\dots,N$ and using (5.1) gives

[TABLE]

and recalling that $W_{i}=S_{i}-\mathbf{i}$ yields

[TABLE]

Since $\tau-\tau^{2}/2$ is clearly maximised at $\tau=1$, the proof is complete. ∎

5.3.1 Proof of Theorem 3

We will prove the theorem by establishing the following facts:

The sequence $\|F^{\prime}(\gamma_{j})\|$ converge to zero as $j\to\infty$. 2. 2.

The sequence $\{\gamma_{j}\}$ is stays in a compact subset of $\mathcal{P}_{2}(\mathbb{R}^{d})$. 3. 3.

The mapping $\gamma\mapsto\|F^{\prime}(\gamma)\|^{2}$ is continuous.

The first two are relatively straightforward, and are proven in the form of the following two Lemmas.

Lemma 3**.**

The objective value of the Fréchet functional decreases at each step of Algorithm 1, and $\|F^{\prime}(\gamma_{j})\|$ vanishes as $j\to\infty$.

Proof.

The first statement is clear from Lemma 2, from which it also follows that

[TABLE]

Consequently, the series at the left-hand side converges whence $\|F^{\prime}(\gamma_{j})\|^{2}\to 0$. ∎

Lemma 4**.**

The sequence generated by Algorithm 1 stays in a compact subset of the Wasserstein space $\mathcal{P}_{2}(\mathbb{R}^{d})$.

Proof.

For any $\epsilon>0$ there exists a compact convex set $K_{\epsilon}$ such that $\mu^{i}(K_{\epsilon})>1-\epsilon/N$ for $i=1,\dots,N$. Let $A^{i}=(\mathbf{t}_{\gamma_{j}}^{\mu^{i}})^{-1}(K_{\epsilon})$, $A=\cap_{i=1}^{N}A^{i}$. Then $\gamma_{j}(A^{i})>1-\epsilon/N$, so that $\gamma_{j}(A)>1-\epsilon$. Since $K_{\epsilon}$ is convex, $T_{j}(x)\in K_{\epsilon}$ for any $x\in A$, so that

[TABLE]

We shall now show that any weakly convergent subsequence of $\{\gamma_{j}\}$ is in fact convergent in the Wasserstein space. By Theorem 7.12 in Villani [74], it suffices to show that

[TABLE]

For simplicity, we shall show this under the stronger assumption that the measures $\mu^{1},\dots,\mu^{N}$ have a finite third moment

[TABLE]

In Section 8 we show that (5.2) holds even if (5.3) does not.

For any $j\geq 1$ it holds that

[TABLE]

This implies that for any $R>0$ and any $j>0$,

[TABLE]

and (5.2) follows. ∎

The third statement (continuity of the gradient) is much more subtle to establish. We will prove it in two steps: first we establish a Proposition, giving sufficient conditions for the third statement to hold true. Then, we will verify that the conditions of the Proposition are satisfied in the setting of Theorem 3.3, in the form of a Lemma and a Corollary. We start with the proposition.

Proposition 3** (Continuity of $F^{\prime}$).**

Let $\mu^{1},\dots,\mu^{N}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be given regular measures, and consider a sequence $\gamma_{n}$ of regular measures that converges in $\mathcal{P}_{2}(\mathbb{R}^{d})$ to a regular measure $\gamma$. If the densities of $\gamma_{n}$ are uniformly bounded, then $\|F^{\prime}(\gamma_{n})\|^{2}\to\|F^{\prime}(\gamma)\|^{2}$.

Proof.

The regularity of $\gamma_{n}$ and $\gamma$ implies that $F$ is indeed differentiable there, and so it needs to be shown that

[TABLE]

Denote the integrands by $g_{n}$ and $g$ respectively. At a given $x\in\mathbb{R}^{d}$, $g_{n}(x)$ can be undefined, either because some $\mathbf{t}_{\gamma_{n}}^{\mu^{i}}(x)$ is empty, or because they can be multivalued. Redefine $g_{n}(x)$ at such points by setting it to 0 in the former case and choosing an arbitrary representative otherwise. Since the set of these ambiguity points is a $\gamma_{n}$-null set (because $\gamma_{n}$ is absolutely continuous), this modification does not affect the value of the integral ${\int\!g_{n}\,\mathrm{d}\gamma_{n}}$. Apply the same procedure to $g$. Then $g_{n}$ and $g$ are finite and nonnegative throughout $\mathbb{R}^{d}$. Absolute continuity of $\gamma$, Remark 2.3 in [3] and Proposition 5 imply together that the set of points where $g$ is not continuous is a $\gamma$-null set.

Next, we approximate $g_{n}$ and $g$ by bounded functions as follows. Since $\gamma_{n}$ converge in the Wasserstein space, they satisfy (5.2) by [74, Theorem 7.12]. It is easy to see that this implies the uniform absolute continuity

[TABLE]

The $\delta$’s can be chosen in such a way that (5.4) holds true for the finite collection $\{\mu^{1},\dots,\mu^{N}\}$ as well. Fix $\epsilon>0$, set $\delta=\delta_{\epsilon}$ as in (5.4), and let $A_{n}=\{x:g_{n}(x)\geq 4R\}$, where $R=R_{\epsilon}\geq 1$ is such that (using (5.2))

[TABLE]

The bound

[TABLE]

implies that

[TABLE]

To deal with the sets in the union observe that (since $\mathbf{t}_{\gamma_{n}}^{\mu^{i}}$ is $\gamma_{n}$-almost surely injective),

[TABLE]

so that $\gamma_{n}(A_{n})<\delta$. We use this in conjunction with (5.4) to bound

[TABLE]

where we have used the measure-preservation property $\mu^{i}(\mathbf{t}_{\gamma_{n}}^{\mu^{i}}(A_{n}))=\gamma_{n}(A_{n})<\delta$.

Define the truncation $g_{n,R}(x)=\min(g_{n}(x),4R)$. Then $0\leq g_{n}-g_{n,R}\leq g_{n}\mathbf{1}\{g_{n}>4R\}$, so

[TABLE]

The analogous truncated function $g_{R}$ satisfies

[TABLE]

Let $E=\mathrm{supp}(\gamma)$. Proposition 6 (Section 5.5) implies pointwise convergence of $\mathbf{t}_{\gamma_{n}}^{\mu^{i}}(x)$ to $\mathbf{t}_{\gamma}^{\mu^{i}}(x)$ for any $i=1,\dots,N$ and any $x\in E\setminus\mathcal{N}$, where $\mathcal{N}=\cup_{i=1}^{N}\mathcal{N}^{i}$ and

[TABLE]

Thus, $g_{n}$ and $g$ are univalued functions defined throughout $\mathbb{R}^{d}$, and $g_{n}\to g$ pointwise on $x\in E\setminus\mathcal{N}$ (for whatever choice of representatives selected to define $g_{n}$); consequently, $g_{n,R}\to g_{R}$ on $E\setminus\mathcal{N}$.

In order to restrict the integrands to a bounded set we invoke the tightness of the sequence $(\gamma_{n})$ and introduce a compact set $K_{\epsilon}$ such that $\gamma_{n}(\mathbb{R}^{d}\setminus K_{\epsilon})<\epsilon/R$ for all $n$. Clearly, $g_{n,R}\to g_{R}$ on $E^{\prime}=K_{\epsilon}\cap E\setminus\mathcal{N}$, and by Egorov’s theorem (valid as $\mathrm{Leb}(E^{\prime})\leq\mathrm{Leb}(K_{\epsilon})<\infty$), there exists a Borel set $\Omega=\Omega_{\epsilon}\subseteq E^{\prime}$ on which the convergence is uniform, and $\mathrm{Leb}(E^{\prime}\setminus\Omega)<\epsilon/R$. Let us write

[TABLE]

and bound each of the three integrals at the right-hand side as $n\to\infty$.

The first integral vanishes as $n\to\infty$, by (5.5) and the Portmanteau lemma (Lemma 9, Section 5.5). For a given $\Omega$, the second integral vanishes as $n\to\infty$, since $g_{n,R}$ converge to $g_{R}$ uniformly. The third integral is bounded by $8R\gamma_{n}(\mathbb{R}^{d}\setminus\Omega)$. The latter set is a subset of $\mathcal{N}\cup(E^{\prime}\setminus\Omega)\cup(\mathbb{R}^{d}\setminus E)\cup(\mathbb{R}^{d}\setminus K_{\epsilon})$, where the first set is Lebesgue-negligible and the second has Lebesgue measure smaller than $\epsilon/R$. The hypothesis of the densities of $\gamma_{n}$ implies that $\gamma_{n}(A)\leq C\mathrm{Leb}(A)$ for any Borel set $A\subseteq\mathbb{R}^{d}$ and any $n\in\mathbb{N}$; it follows from this and $\gamma_{n}(\mathbb{R}^{d}\setminus K_{\epsilon})<\epsilon/R$ that

[TABLE]

Write the open set $E_{1}=\mathbb{R}^{d}\setminus E$ as a countable union of closed sets $A_{k}$ with $\mathrm{Leb}(E_{1}\setminus A_{k})<1/k$, and conclude that

[TABLE]

where we have used the Portmanteau lemma again, $A_{k}\cap\mathrm{supp}(\gamma)=\emptyset$ and $\gamma_{n}(A)\leq C\mathrm{Leb}(A)$. Consequently, for all $k$

[TABLE]

Letting $k\to\infty$, then incorporating the truncation error yields

[TABLE]

The proof is complete upon noticing that $\epsilon$ is arbitrary. ∎

Our proof will now be complete if we show that the sequence $\gamma_{k}$ generated by the algorithm satisfies the assumptions of the last Proposition. First we show that limits of the sequence are indeed regular.

Proposition 4** (Sequence has bounded density).**

Let $\mu^{i}$ have density $g^{i}$ for $i=1,\dots,N$ and let $\gamma_{0}$ be a regular probability measure. Then the density of $\gamma_{1}$ is bounded by a constant $C_{\mu}=\min\{N^{d-1}\max_{i}\|g^{i}\|_{\infty},N^{d}\min_{i}\|g^{i}\|_{\infty}\}$ that depends only on $\{\mu^{1},\dots,\mu^{N}\}$.

Proof.

Let $h_{i}$ be the density of $\gamma_{i}$. By the change of variables formula, for $\gamma_{0}$-almost any $x$

[TABLE]

Fiedler [31] shows that if $B_{1}$ and $B_{2}$ are $d\times d$ positive semidefinite matrices with eigenvalues $0\leq\alpha_{i},\beta_{i}$, then

[TABLE]

The right-hand side contains $2^{d}$ nonnegative summands of which two are $\mathrm{det}B_{1}$ and $\mathrm{det}B_{2}$, and so we see that $\mathrm{det}(B_{1}+B_{2})\geq\mathrm{det}B_{1}+\mathrm{det}B_{2}$. (One can show the stronger result $\sqrt[d]{\mathrm{det}(B_{1}+B_{2})}\geq\sqrt[d]{\mathrm{det}B_{1}}+\sqrt[d]{\mathrm{det}B_{2}}$.) Since $\nabla\mathbf{t}_{\gamma_{0}}^{\gamma_{1}}$ is an average of $N$ $d\times d$ positive semidefinite matrices, we obtain

[TABLE]

Let $\Sigma$ be the set of points where this inequality holds; then $\gamma_{0}(\Sigma)=1$. Hence

[TABLE]

Thus $\gamma_{1}$-almost surely,

[TABLE]

For $C_{\mu}$ to be finite it suffices that $\|g^{i}\|_{\infty}$ be finite for some $i$. ∎

Our task is now essentially complete. All that remains is to show:

Corollary 5** (Limits are regular).**

Every limit of the sequence generated by Algorithm 1 is absolutely continuous provided the density of $\mu^{i}$ is bounded for some $i$.

Proof.

Each $\gamma_{k}$ ($k=1,2,\dots$) has a density that is bounded by the finite constant $C_{\mu}$. For any open set $O$, $\liminf\gamma_{k}(O)\leq C_{\mu}\mathrm{Leb}(O)$, so any limit point $\gamma$ of $(\gamma_{k})$ is such that $\gamma(O)\leq C_{\mu}\mathrm{Leb}(O)$ by the Portmanteau lemma. It follows that $\gamma$ is absolutely continuous with density bounded by $C_{\mu}$. We note that Agueh and Carlier [2] show that the density of the Fréchet mean is bounded by $N^{d}\min_{i}\|g^{i}\|_{\infty}\geq C_{\mu}$, a slightly weaker bound. ∎

Proof of Theorem 4.

Let $E=\mathrm{supp}(\bar{\mu})$ and set $A^{i}=E^{\mathrm{den}}\cap\{x:\mathbf{t}_{\bar{\mu}}^{\mu^{i}}(x)\textrm{ is multivalued}\}$. By Corollary 6 $\bar{\mu}(A^{i})=1$. Choose $A=\cap_{i=1}^{N}A^{i}$ and apply Proposition 6. This proves the first assertion.

Now let $E^{i}=\mathrm{supp}(\mu^{i})$ and set $B^{i}=(E^{i})^{\mathrm{den}}\cap\{x:\mathbf{t}_{\mu^{i}}^{\bar{\mu}}(x)\textrm{ is univalued}\}$. Since $\mu^{i}$ is regular, $\mu^{i}(B^{i})=1$. Apply Proposition 6. If in addition $E^{1}=\dots=E^{N}$ then $\mu^{i}(B)=1$ for $B=\cap B^{i}$. ∎

Proof of Corollary 3.

The proof is very similar to that of Proposition 3. Define $\eta_{j},\eta\in\mathcal{P}_{2}((\mathbb{R}^{d})^{N+1})$ by

[TABLE]

We establish convergence of $\eta_{j}$ to $\eta$. Since the optimal multicouplings are marginals of $\eta_{j}$ and $\eta$ their convergence follow. Let $h:(\mathbb{R}^{d})^{N+1}\to\mathbb{R}$ be any continuous function such that

[TABLE]

Define $g_{j}:\mathbb{R}^{d}\to\mathbb{R}$ by $g_{j}(x)=h(\mathbf{t}_{\gamma_{j}}^{\mu^{1}},\dots\mathbf{t}_{\gamma_{j}}^{\mu^{n}},x)$ and analogously define $g$. By [74, Theorem 7.12] it suffices to show that (if this holds for $h$, it also holds for $a+bh$ with $a,b$ scalars)

[TABLE]

(In Proposition 3 we had $h=\|y-\bar{t}\|^{2}$.) Since $h$ is continuous, we can modify $g_{n}$ and $g$ to be well-defined, finite and so that $g$ be continuous $\gamma$-almost surely. Define $R$ as in the proof of Proposition 3, $A_{j}=\{x:|g_{j}(x)|\geq 4R\}$, invoke (5.4) and translate the bound on $h$ to a bound on $|g_{j}|$ to conclude that ${\int_{A_{j}}\!|g_{j}|\,\mathrm{d}\gamma_{j}}\leq 4\epsilon$. Carry out the same (now two-sided) truncation $g_{j,R}(x)=\max(-4R,\min(g_{j}(x),4R))$ to obtain $|g_{j}-g_{j,R}|\leq|g_{j}|\mathbf{1}\{|g_{j}|>4R\}$, $|g_{R}|\leq 4R$ and $g_{R}$ is continuous $\gamma$-almost surely (see (5.5)). The rest can be done as in the proof of Proposition 3, since it did not depend on the precise form of $g$. ∎

5.4 Proofs of Statements in Section 4.1

Proof of Theorem 5.

Since $\lambda$ is regular and $T$ is injective with nonsingular derivative, $\Lambda=T\#\lambda$ is also regular by Lemma 5.5.3 in [7]. Moreover, $\Lambda$ is supported on $K$ because $T$ takes values there. Consequentely, the Fréchet mean of $\Lambda$ is unique and supported itself on $K$; this is essentially a consequence of Corollary 2.9 in [5]. For tidiness, we provide the full details in Section 8.

In view of the preceding paragraph, it suffices to show that

[TABLE]

As a gradient of a convex function, $T=\mathbf{t}_{\lambda}^{\Lambda}$ is optimal. Let $\phi$ be the convex potential of $T$, and define $\phi^{*}$ its Legendre transform. Then the pair $(\|x\|^{2}/2-\phi,\|y\|^{2}/2-\phi^{*})$ is dual optimal. Invoking strong duality for $\lambda$ and weak duality for $\theta$, we find

[TABLE]

By Fubini’s theorem (see the supplementary material for a justification), we have

[TABLE]

The function $\mathbb{E}T$ is continuous (by the bounded convergence theorem and boundedness of $K$), so equals the identity for all $x\in K$. Again by Fubini’s theorem (see the supplementary material), it follows that $\mathbb{E}\phi(x)=\|x\|^{2}/2$ for all $x\in K$, perhaps up to an additive constant. Since $\theta$ and $\lambda$ are both supported on $K$, the integrals with respect to $\lambda$ and $\theta$ vanish, and this completes the proof. ∎

As part of our proofs, we will need to control the Wasserstein distance between the regularised measures and their true counterparts:

Lemma 5**.**

The smooth measure $\widehat{\Lambda}_{i}$ defined by (4.1) satisfies

[TABLE]

where $C_{\psi,K}$ is a (finite) constant that depends only on $\psi$ and $K$.

We prove the lemma in the supplementary material (Section 8).

Remark 5**.**

There is no need for $\psi$ to be isotropic: it is sufficient that merely

[TABLE]

which is satisfied as long as $\psi$ is continuous and strictly positive.

We now remark that a trivial extension of [61, Lemma 3] yields:

Lemma 6** (Number of points per process is $O(\tau_{n})$).**

If $\tau_{n}/\log n\to\infty$, then there exists a constant $C_{\Pi}>0$, depending only on the distribution of the $\Pi$’s, such that

[TABLE]

In particular, there are no empty point processes, so the normalisation is well-defined.

Proof of Theorem 6.

The proof is very similar to that of Theorem 1 in Panaretos & Zemel [61], and we give the details in the supplementary material (Section 8). ∎

Proof of Theorem 7.

The argument is considerably different than the case $d=1$ considered in [61], and brings into play the geometry of convex functions in $\mathbb{R}^{d}$. Let $i$ be a fixed integer and for $n\geq i$ set

[TABLE]

We wish to show that $u_{n}\to u$ uniformly on compact sets, using our knowledge that

[TABLE]

This follows from Proposition 6 below. To verify the conditions, notice that all the measures are supported on $K=E$, a compact and convex set. Furthermore $\mu_{n}$, $\mu$ and $\nu$ all have strictly positive densities there, so their support is exactly $K$. Continuity of $u$ on $\mathrm{int}(K)$ follows from the assumptions that $T_{i}$ and $T_{i}^{-1}$ are continuous. The finiteness in (5.7) follows from the compactness of $K$, and the uniqueness follows from the regularity of $\mu$.

The same proposition can be applied to show convergence of $\widehat{T}_{i}$ to $T_{i}$ uniformly on $\Omega\subseteq\mathrm{int}(K)$: one needs to reverse the roles of $\mu_{n}$ and $\nu_{n}$ and of $\mu$ to $\nu$, and notice that $\nu$ too is regular, which guarantees the uniqueness in (5.7). ∎

Proof of Corollary 4.

The square of the distance is

[TABLE]

and this is well-defined (that is, $\Pi_{i}(K)>0$) almost surely for $n$ large enough by Lemma 6. Since $\lambda(\partial K)=0$, almost surely there are no points on the boundary and the integral can be taken on the interior of $K$. Let $\Omega\subseteq\mathrm{int}(K)$ be compact and split the integral to $\Omega$ and its complement. Then

[TABLE]

by the law of large numbers. Since the interior of $K$ can be written as a countable union of compact sets, the right-hand side can be made arbitrarily small by selection of $\Omega$.

Let us now consider the integral on $\Omega$. Since

[TABLE]

and $T_{i}(\Omega)$ is compact, we only need to show that it is included in $\mathrm{int}(K)$ in order to apply Theorem 7. Suppose towards contradiction that $y=T_{i}(x)\in\partial K$ for $x\in\mathrm{int}(K)$. Let $\alpha\in\mathbb{R}^{d}\setminus\{0\}$ with $\langle y,\alpha\rangle\geq\sup\langle K,\alpha\rangle$. Let $x^{\prime}=x+t\alpha$ for $t>0$ small enough such that $x^{\prime}\in\mathrm{int}(K)$. Then $y^{\prime}=T_{i}(x^{\prime})\in K$, so that

[TABLE]

Either condition in the statement of the corollary imply that $y^{\prime}=y$, in contradiction to $T_{i}$ being injective. ∎

5.5 Monotone Operators, Optimal Transportation, Stochastic Convergence

This section contains the statements and proofs of analytical results needed in our proofs, culminating in Proposition 6. The latter is the backbone result needed for the proofs of Theorem 7, Theorem 3 (more precisely, Proposition 3) and Theorem 4. Rather than start with all the background definitions we will define the necessary objects en route.

We shall follow the notation and terminology of Alberti & Ambrosio [3]. Let $u$ be a set-valued function (or multifunction) on $\mathbb{R}^{d}$, that is, $u:\mathbb{R}^{d}\to 2^{\mathbb{R}^{d}}$. It is said that $u$ is monotone if

[TABLE]

When $d=1$, the definition reduces to $u$ being a nondecreasing (set-valued) function. It is said that $u$ is maximal if no points can be added to its graph while preserving monotonicity:

[TABLE]

We sometimes use the notation $(x,y)\in u$ to mean $y\in u(x)$. Note that $u(x)$ can be empty, even when $u$ is maximal.

The relevance of monotonicity stems from the fact that subdifferentials of convex functions are monotone. That is, if $\varphi:\mathbb{R}^{d}\to\mathbb{R}\cup\{\infty\}$ is lower semicontinuous and convex (and not identically infinite), then $u=\partial\varphi$ is maximally monotone [3, Section 7], where

[TABLE]

is the subdifferential of $\varphi$ at $x$. Here $u(x)=\emptyset$ if $\varphi(x)=\infty$.

We will use extensively the continuity of $u$ at points where it is univalued.

Proposition 5** (Continuity at Singletons).**

Let $u$ be a maximal monotone function, and suppose that $u(x)=\{y\}$ is a singleton. Then $u$ is nonempty on some neighbourhood of $x$ and it is continuous at $x$: if $x_{n}\to x$ and $y_{n}\in u(x_{n})$, then $y_{n}\to y$.

Proof.

See [3, Corollary 1.3(4)]. Notice that this result implies that differentiable convex functions are continuously differentiable [66, Corollary 25.5.1].∎

It turns out that when $u$ is univalued, monotonicity is a local property. To state the result in the general form that we shall use, we need to introduce the notion of points of Lebesgue density.

Let $B_{r}(y)=\{x:\|x-y\|<r\}$ for $r\geq 0$ and $y\in\mathbb{R}^{d}$. A point $x_{0}$ is of Lebesgue density of a measurable set $G\subseteq\mathbb{R}^{d}$ if for any $\epsilon>0$ there exists $t_{\epsilon}>0$ such that

[TABLE]

We denote the set of points of Lebesgue density of $G$ by $G^{\rm{den}}$. Clearly, $G^{\rm{den}}$ lies between $\mathrm{int}(G)$ and $\overline{G}$. Stein and Shakarchi [72, Chapter 3, Corollary 1.5] show that almost any point of $G$ is in $G^{\rm{den}}$. By the Hahn–Banach theorem, $G^{\rm{den}}\subseteq\mathrm{int}(\mathrm{conv}(G))$.

Lemma 7** (Density Points and Distance).**

Let $x_{0}$ be a point of Lebesgue density of a measurable set $G\subseteq\mathbb{R}^{d}$. Then

[TABLE]

This result was given as an exercise in [72]; for completeness we provide a full proof in the supplementary material (Section 8).

Lemma 8** (Local Monotonicity).**

Let $u$ be a maximal monotone function such that $u(x_{0})=\{y_{0}\}$. Suppose that $x_{0}$ is a point of Lebesgue density of a set $G$ satisfying

[TABLE]

Then $y^{*}=y_{0}$. In particular, the result is true if the inequality holds on $G=O\setminus\mathcal{N}$ with $\emptyset\neq O$ open and $\mathcal{N}$ Lebesgue negligible.

Proof.

Set $z_{t}=x_{0}+t(y^{*}-y_{0})$ for $t>0$ small. It may be that $z_{t}\notin G$; but Lemma 7 guarantees existence of $x_{t}\in G$ with $\|x_{t}-z_{t}\|/t\to 0$. By Proposition 5 $u(x_{t})$ is nonempty for $t$ small enough. For $y_{t}\in u(x_{t})$,

[TABLE]

Rearrangement, division by $t>0$ and application of the Cauchy–Schwartz inequality gives

[TABLE]

As $t\searrow 0$ the right-hand side vanishes, since $y_{t}\to y_{0}$ (Proposition 5) and $\|x_{t}-z_{t}\|/t\to 0$. It follows that $y^{*}=y_{0}$. ∎

This concludes the necessary discussion on monotone operators. We will now state some necessary results on optimal transportation maps, and specifically their convergence properties. Consider the following setting: let $\{\mu_{n}\}$, $\{\nu_{n}\}$ be two sequences of probability measures on $\mathbb{R}^{d}$ that converge weakly to $\mu$ and $\nu$ respectively. Let $\pi_{n}$ be an optimal coupling between $\mu_{n}$ and $\nu_{n}$ having finite cost, which is supported on the graph of a subdifferential of a proper (not identically infinite) convex lower semicontinuous function $\varphi_{n}$ [74, Chapter 2]. The set-valued function $u_{n}=\partial\varphi_{n}$ that maps $x$ to the subdifferential of $\varphi_{n}$ at $x$ is maximally monotone [3, Section 7]. The appropriate functions for $\mu$ and $\nu$ will be denoted by $\varphi$ and $u=\partial\varphi$ and the optimal coupling by $\pi$. This setting will be succinctly referred to by the equation

[TABLE]

We notice now that uniqueness of $\pi$ and the stability of optimal transportation imply that $\pi_{n}$ converge weakly to $\pi$ (even if $\pi_{n}$ is not unique); see Schachermayer & Teichmann [69, Theorem 3] or Cuesta-Albertos et al. [25, Theorem 3.2]. This weak convergence will be used in the following form:

Lemma 9** (Portmanteau).**

Weak convergence of Borel probability measures $\mu_{k}$ to $\mu$ on $\mathbb{R}^{d}$ is equivalent to any of the following conditions:

• (I)

for any open set $G$, $\liminf\mu_{k}(G)\geq\mu(G)$;

• (II)

for any closed set $F$, $\limsup\mu_{k}(F)\leq\mu(F)$;

• (III)

${\int\!h\,\mathrm{d}\mu_{k}}\to{\int\!h\,\mathrm{d}\mu}$* for any bounded measurable $h$ whose set of discontinuity points is a $\mu$-null set.*

Proof.

The equivalence with the first two conditions is classical and can be found in Billingsley [15, Theorem 2.1]; for the third, see Pollard [64, Section III.2]. ∎

We shall now translate this into convergence of $u_{n}$ to $u$ under certain regularity conditions.

Proposition 6** (Uniform Convergence of Optimal Maps).**

In the setting of Display (5.7), denote $E=\mathrm{supp}(\mu)$.

Let $\Omega$ be a compact subset of $E^{\mathrm{den}}$ on which $u$ is univalued, where $E^{\mathrm{den}}$ is the set of points of Lebesgue density of $E$. Then $u_{n}$ converges to $u$ uniformly on $\Omega$: $u_{n}(x)$ is nonempty for all $x\in\Omega$ and all $n>N_{\Omega}$, and

[TABLE]

In particular, if $u$ is univalued throughout $\mathrm{int}(E)$ (so that $\varphi\in C^{1}$ there), then uniform convergence holds for any compact $\Omega\subset\mathrm{int}(E)$.

Corollary 6** (Pointwise convergence $\mu$-almost surely).**

If in addition $\mu$ is absolutely continuous then $u_{n}(x)\to u(x)$ $\mu$-almost surely.

Proof.

The set of points $x\in E$ for which $\Omega=\{x\}$ fails to satisfy the conditions of Proposition 6 is included in

[TABLE]

(Since $u$ is nonempty on $\mathrm{int}(\mathrm{conv}(E))$ by [3, Corollary 1.3(2)].) Both sets are Lebesgue-negligible (see [3, Remark 2.3] for the latter), and $\mu$ is absolutely continuous. ∎

Remark 6**.**

In the setting of Theorem 7, $E$ is convex, $\mu$ is absolutely continuous, and $u$ is univalued on $\mathrm{int}(E)$, so one can take any $\Omega\subseteq\mathrm{int}(E)$, without the need to introduce Lebesgue density. The more general statement of the proposition is used in the proof of Proposition 3, where we have no control on the support of $\gamma$ or the regularity of the transport maps.

We split the proof of Proposition 6 into two steps: (1) Limit points of the graphs of $u_{n}$ are in the graph of $u$ (Lemma 11); (2) Points in the graphs of $u_{n}$ stay in a bounded set (Proposition 7). Each of these points will be proven using one intermediate lemma.

Lemma 10** (Points in the limit graph are limit points).**

Assume (5.7). For any $x_{0}\in\mathrm{supp}(\mu)$ such that $u(x_{0})=\{y_{0}\}$ is a singleton there exists a subsequence $(x_{n_{k}},y_{n_{k}})\in u_{n_{k}}$ that converges to $(x_{0},y_{0})$.

Proof.

Since $u=\partial\varphi$ is a maximal monotone function [3, Section 7] that is univalued at $x_{0}$, it is continuous there (Proposition 5). This means that for any $\epsilon>0$ there exists $\delta>0$ such that if $x\in B_{\delta}(x_{0})=\{x:\|x-x_{0}\|<\delta\}$ then $u(x)$ is nonempty and if $y\in u(x)$, then $\|y-y_{0}\|<\epsilon$. Take $\epsilon_{k}\to 0$ and corresponding $\delta_{k}\to 0$, and set $B_{k}=B_{\delta_{k}}(x_{0})$, $V_{k}=B_{\epsilon_{k}}(y_{0})$. Then $u(B_{k})\subseteq V_{k}$, so

[TABLE]

because $B_{k}$ is a neighbourhood of $x_{0}\in\mathrm{supp}(\mu)$. Since $B_{k}\times V_{k}$ is open, we have by the Portmanteau lemma that $\pi_{n}(B_{k}\times V_{k})>0$ for $n$ large. Consequently, there exists $n_{k}$ such that

[TABLE]

Since $\pi_{n_{k}}$ is concentrated on the graph of $u_{n_{k}}$, it follows that there exist $(x_{n_{k}},y_{n_{k}})\in u_{n_{k}}$ with $\|x_{n_{k}}-x_{0}\|<\delta_{k}$ and $\|y_{n_{k}}-y_{0}\|<\epsilon_{k}$. Hence $(x_{n_{k}},y_{n_{k}})\to(x_{0},y_{0})$. ∎

Lemma 11** (Limit points are in the limit graph).**

Assume that (5.7) holds and denote $E=\mathrm{supp}(\mu)$. If a subsequence $(x_{n_{k}},y_{n_{k}})\in u_{n_{k}}$ converges to $(x_{0},y^{*})$, where $x_{0}$ is a point of Lebesgue density of $E$, and $u(x_{0})$ is a singleton, then $y^{*}=u(x_{0})$. In particular, the statement is true if $x_{0}\in\mathrm{int}(E)$ and $u(x_{0})$ is a singleton.

Proof.

The set $\mathcal{N}\subseteq\mathbb{R}^{d}$ of points where $u$ contains more than one element is Lebesgue negligible [3, Remark 2.3]. There exists a neighbourhood $V$ of $x_{0}$ on which $u$ is nonempty (Proposition 5). Thus, $x_{0}$ is a point of Lebesgue density of $G=(E\cap V)\setminus\mathcal{N}$, and $u(x)$ is a singleton for every $x\in G$. Fix such an $x$ and set $y=u(x)$. By Lemma 10 (applied to $\{u_{n_{k}}\}_{k=1}^{\infty}$ at $x$) there exist sequences $x_{n_{k_{l}}}^{\prime}\to x$ and $y_{n_{k_{l}}}^{\prime}\to y$ with $(x_{n_{k_{l}}}^{\prime},y_{n_{k_{l}}}^{\prime})\in u_{n_{k_{l}}}$. Consequently,

[TABLE]

This holds for any $(x,y)\in u$ such that $x\in G$. Since $x_{0}$ is a point of Lebesgue density of $G$ (and $u$ is maximal), it follows from Lemma 8 that $y^{*}=u(x_{0})$. ∎

Let $B_{\epsilon}^{\infty}(x_{0})=\{x:\|x-x_{0}\|_{\infty}<\epsilon\}$ be the $\ell_{\infty}$ ball around $x_{0}$ and $\overline{B}_{\epsilon}^{\infty}(x_{0})$ its closure.

Lemma 12** (Continuity of Convex Hulls).**

Let $Z=\{z_{i}\}\subseteq\mathbb{R}^{d}$ be a set of points whose convex hull, $\mathrm{conv}(Z)$, includes $B_{\rho}^{\infty}(x_{0})$ and let $\tilde{Z}=\{\tilde{z}_{i}\}$ be a set of points such that $\|\tilde{z}_{i}-z_{i}\|_{\infty}\leq\epsilon$. Then the convex hull of $\tilde{Z}$ includes $B_{\rho-\epsilon}^{\infty}(x_{0})$.

For a proof, see Section 8.

Proposition 7** (Boundedness).**

Suppose that (5.7) holds, and fix a compact $\Omega\subseteq\mathrm{int}(\mathrm{conv}(\mathrm{supp}(\mu)))$. Then for $n>N(\Omega)$ sufficiently large, $u_{n}(x)$ is nonempty for all $x\in\Omega$ and $u_{n}(\Omega)$ is bounded uniformly.

Proof.

Denote $E=\mathrm{supp}(\mu)$ and its convex hull by $F=\mathrm{conv}(E)$. There exists $\delta=\delta(\Omega)>0$ such that the closed $\ell_{\infty}$-ball, $\overline{B}^{\infty}_{3\delta}(\Omega)$, is included in $\mathrm{int}(F)$. Cover $\Omega$ by a finite union of $B_{\delta}^{\infty}(\omega_{j})$, and denote by $Q$ be the finite set of vertices of $\cup_{j}\overline{B}^{\infty}_{3\delta}(\omega_{j})$. Since $Q$ is included in the convex hull of $E$, each point in $Q$ can be written as a convex combination of elements of $E$. We conclude that there exists a finite set $Z=\{z_{1},\dots,z_{m}\}$ of points in $E$ whose convex hull includes $B_{3\delta}^{\infty}(\omega_{j})$ for any $j$.

Let $B_{i}=B_{\delta}^{\infty}(z_{i})$. Since $B_{i}$ is an open neighbourhood of $z_{i}\in E=\mathrm{supp}(\mu)$, the Portmanteau lemma implies that when $n$ is large, $\mu_{n}(B_{i})>\epsilon_{i}=\mu(B_{i})/2$ for any $i=1,\dots,m$. Let $\epsilon=\min_{i}\epsilon_{i}>0$. Since $\{\nu_{n}\}$ is a tight sequence, there exists a compact set $K_{\epsilon}$ such that $\nu_{n}(K_{\epsilon})>1-\epsilon$ for any integer $n$. In particular, there exist $x_{ni}\in B_{i}$ and $y_{ni}\in u_{n}(x_{ni})$ such that $y_{ni}\in K_{\epsilon}$. Application of Lemma 12 to

[TABLE]

and noticing that by definition $\|x_{ni}-z_{i}\|_{\infty}\leq\delta$ yields

[TABLE]

For each $\omega\in\Omega$ there exists $j$ such that $\|\omega-\omega_{j}\|_{\infty}\leq\delta$, so that $\mathrm{conv}(X_{n})\supseteq B_{\delta}^{\infty}(\omega)\supseteq B_{\delta}(\omega)$, since $\ell_{2}$-balls are smaller than $\ell_{\infty}$-balls. Summarising: $\mathrm{conv}(X_{n})\supseteq B_{\delta}(\Omega)$.

By [3, Lemma 1.2(4)] it follows that for any $\omega\in\Omega$ and any $y_{0}\in u_{n}(\omega)$,

[TABLE]

Now observe that the infimum at the right-hand side is bounded by $\|y_{ni}\|\leq\sup_{y\in K_{\epsilon}}\|y\|$. Furthermore, $\|x_{nk}-x_{nl}\|\leq 2\sqrt{d}\delta+\|z_{k}-z_{l}\|$. Hence

[TABLE]

and the right-hand side is independent of $n$. We may therefore conclude that for $n$ large enough, $u_{n}(\Omega$) stays in a compact set; it is nonempty by [3, Corollary 1.3(2)]. ∎

Proof of Proposition 6.

By Proposition 7 when $n>N_{\Omega}$ is large, $u_{n}(x)\neq\emptyset$ for all $x\in\Omega$ and

[TABLE]

where $C_{\Omega,d}$ is a constant that depends only on $\Omega$ (and the dimension $d$).

Suppose that the converse is true, and uniform convergence does not hold. Then there exist $\epsilon>0$ and subsequences $y_{n_{k}}\in u_{n_{k}}(x_{n_{k}})$ such that $x_{n_{k}}\in\Omega$ and

[TABLE]

The $x_{n_{k}}$’s lie in the compact set $\Omega$, whereas by Proposition 7 the $y_{n_{k}}$’s lie in the ball of radius $C_{\Omega,d}$ centred at the origin. Therefore, up to the extraction of a subsequence, we have $x_{n_{k}}\to x\in\Omega$ and $y_{n_{k}}\to y$. By Lemma 11, $y=u(x)$. But $u$ is continuous at $x$ (Proposition 5), whence

[TABLE]

a contradiction. ∎

6 Some Examples

As an illustration, we implement Algorithm 1 in several settings for which pairwise optimal maps can be calculated explicitly at every iteration, allowing for fast computation without error propagation. Indeed, these settings allow for stronger convergence statements to be made on a case-by-case basis. More details on the calculations and properties of each individual scenario are given in Section 8.

6.1 The case $d=1$

When the measures are supported on the real line, the optimal maps have the explicit expression given in Equation (2.1) and one may apply Algorithm 1 starting from one of these measures. Figure 1 plots $N=4$ univariate densities and the Fréchet mean yielded by the algorithm in two different scenarios. At the left, the densities were generated as

[TABLE]

with $\phi$ the standard normal density, and the parameters generated independently as

[TABLE]

At the right of Figure 1, we used a mixture of a shifted gamma and a Gaussian:

[TABLE]

with

[TABLE]

The resulting Fréchet mean density for both settings is shown in thick light blue, and can be seen to capture the bimodal nature of the data. Even though the Fréchet mean of Gaussian mixtures is not a Gaussian mixture itself, it is approximately so, provided that the peaks are separated enough. Figure 8(a) shows the Procrustes maps pushing the Fréchet mean $\bar{\mu}$ to the measures $\mu^{1},\dots,\mu^{N}$ in each case. If one ignores the “middle part” of the $x$ axis, the maps appear (approximately) affine for small values of $x$ and for large values of $x$, indicating how the peaks are shifted. In the middle region, the maps need to “bridge the gap” between the different slopes and intercepts of these affine maps.

6.2 Independence

We next take measures $\mu^{i}$ on $\mathbb{R}^{2}$, having independent marginal densities $f_{X}^{i}$ as in (6.1), and $f_{Y}^{i}$ as in (6.2). Figure 2 shows the density plot of $N=4$ such measures, constructed as the product of the measures from Figure 1. One can distinguish the independence by the “parallel” structure of the figures: for every pair $(y_{1},y_{2})$, the ratio $g(x,y_{1})/g(x,y_{2})$ does not depend on $x$ (and vice versa, interchanging $x$ and $y$). Figure 3 plots the density of the resulting Fréchet mean. We observe that the Fréchet mean captures the four peaks, and their location. Furthermore, the parallel nature of the figure is preserved in the Fréchet mean. Indeed, we prove in the supplement (Section 8) that, unsurprisingly, the Fréchet mean is a product measure.

6.3 Common Copulas

Let $\mu^{i}$ be a measure on $\mathbb{R}^{2}$ with density

[TABLE]

where $f_{X}^{i}$ and $f_{Y}^{i}$ are random densities on the real line with distribution functions $F_{X}^{i}$ and $F_{Y}^{i}$, and $c$ is a copula density. Figure 4 shows the density plot of $N=4$ such measures, with $f_{X}^{i}$ generated as in (6.1), $f_{Y}^{i}$ as in (6.2), and $c$ is the Frank($-8$) copula density, while Figure 5 plots the density of the Fréchet mean obtained. (For ease of comparison we use the same realisations of the densities that appear in Figure 1.) The Fréchet mean can be seen to preserve the shape of the density, having four clearly distinguished peaks. Figure 8(b), depicting the resulting Procrustes maps, allows for a clearer interpretation: for instance the leftmost plot (in black) shows more clearly that the map splits the mass around $x=-2$ to a much wider interval; and conversely a very large amount mass is sent to $x\approx 2$. This rather extreme behaviour matches the peak of the density of $\mu^{1}$ located at $x=2$.

The first three scenarios are examples of situations where the measures $\{\mu^{i}\}$ are compatible with each other in the sense that $\mathbf{t}^{\mu^{k}}_{\mu^{j}}\circ\mathbf{t}_{\mu^{i}}^{\mu^{j}}=\mathbf{t}_{\mu^{i}}^{\mu^{k}}$. Boissard et al. [16] tackle the problem of finding the Fréchet mean in such a setting, by means of the iterated barycentre. In the supplementary material (Section 8) we show that Algorithm 1 will always converges to the Fréchet mean, provided the initial point $\gamma_{0}$ is compatible with $\{\mu^{i}\}$ (for instance, if $\gamma_{0}=\mu^{i}$). In fact, we show that convergence is established after a single iteration of the algorithm. Since optimal maps are gradients of convex potentials, they must have positive definite derivatives. Under regularity conditions, compatibility is essentially equivalent to the commutativity of the $d\times d$ matrices $\nabla\mathbf{t}_{\mu^{j}}^{\mu^{k}}(\mathbf{t}_{\mu^{i}}^{\mu^{j}}(x))$ and $\nabla\mathbf{t}_{\mu^{i}}^{\mu^{j}}(x)$ for $\mu^{i}$-almost any $x$. We next discuss examples where this condition fails.

6.4 Gaussian measures

Suppose that each $\mu^{i}$ follows a non-degenerate multivariate Gaussian distribution with mean [math] and covariance matrix $S_{i}$. The optimal maps are known to be linear and admit the explicit formula (Dowson & Landau [29]; Olkin & Pukelsheim [60])

[TABLE]

If the initial point $\gamma_{0}$ is another Gaussian measure with covariance matrix $\Gamma_{0}$, then by the linearity of the maps one sees that $\gamma_{k}\sim\mathcal{N}(0,\Gamma_{k})$ for some positive definite $\Gamma_{k}$. Thus, one can calculate the optimal maps at each iteration; in the supplement (Section 8) we prove that $\gamma_{k}$ must converge to the unique Fréchet mean, which is also a Gaussian measure. This example is also studied independently in Álvarez-Esteban et al. [6, Section 4], where an alternative proof can be found. Our proof is shorter and arguably simpler, but the proof in [6] shows the additional property that the traces of the matrix iterates are monotonically increasing.

Notice that the Gaussian measures $\{\mu^{i}\}$ will be compatible if $S_{i}S_{j}=S_{j}S_{i}$, but they might well fail to be. Thus, the algorithm does not converge in one step. We observed, however, rapid convergence of the iterates of Algorithm 1 to the Fréchet mean, even for rather large values of $N$ and $d$. Figure 6 shows density plots of $N=4$ centred Gaussian measures on $\mathbb{R}^{2}$ with covariances $S_{i}\sim\mathrm{Wishart}(I_{2},2)$, and Figure 7 shows the density of the resulting Fréchet mean. In this particular example, the algorithm needed 11 iterations starting from the identity matrix. The corresponding Procrustes registration maps are displayed in Figure 8(c). It is apparent from the figure that these maps are linear, and after a more careful reflection one can be convinced that their average is the identity. The four plots in the figure are remarkably different, in accordance with the measures themselves having widely varying condition numbers and orientations; $\mu^{3}$ and more so $\mu^{4}$ are very concentrated, so the registration maps “sweep” the mass towards zero. In contrast, the registration maps to $\mu^{1}$ and $\mu^{2}$ spread the mass out away from the origin.

6.5 Partially Gaussian Trivariate Measures

We now apply Algorithm 1 in a situation that entangles two of the previous settings. Let $U$ be a $3\times 3$ real orthogonal matrix with columns $U_{1}$, $U_{2}$, $U_{3}$ and let $\mu^{i}$ have density

[TABLE]

with $f^{i}$ bounded density on the real line and $S^{i}\in\mathbb{R}^{2\times 2}$ positive definite. We simulated $N=4$ such densities with $f^{i}$ as in (6.1) and $S^{i}\sim\mathrm{Wishart}(I_{2},2)$. We apply Algorithm 1 to this collection of measures and find their Fréchet mean (in Section 8 we provide precise details on how the optimal maps were calculated). Figure 9 shows level set of the resulting densities for some specific values. The bimodal nature of $f^{i}$ implies that for most values of $a$, $\{x:f^{i}(x)=a\}$ has four elements. Hence the level sets in the figures are unions of four separate parts, with each peak of $f^{i}$ contributing two parts that form together the boundary of an ellipsoid in $\mathbb{R}^{3}$ (see Figure 10). The principal axes of these ellipsoids and their position in $\mathbb{R}^{3}$ differ between the measures, but the Fréchet mean can be viewed as an average of those in some sense.

In terms of orientation (principal axes) of the ellipsoids, the Fréchet mean is most similar to $\mu^{1}$ and $\mu^{2}$, whose orientations are similar to one another.

In the most general examples, one might not be able to analytically obtain the optimal maps at each iteration. In such situations, one needs to resort to numerical schemes such as Benamou & Brenier [10], Haber et al. [40] or Chartrand et al. [23] to obtain the $N$ optimal maps at each iteration (see the concluding remarks for further discussion about numerical issues). Usually such schemes are iterative themselves, so one must take care in managing propagation of errors resulting from using approximate rather than exact transport maps.

7 Concluding Remarks

While the algorithm and the convergence analysis in this work were discussed in the context of absolutely continuous measures, it is worth mentioning the possibility of applying it to discrete measures in some special cases. Specifically, suppose that each measure $\mu^{i}$ is uniform on a set of $M$ distinct points, $\{x^{i}_{m}\}_{m=1}^{M}$. Define as in Anderes et al. [9] the set

[TABLE]

of averages of choices of points from the supports of $\{\mu^{i}\}$. Let $\gamma_{0}$ be an initial measure, uniform on $M$ distinct points as well. There exist optimal maps (not necessarily unique) from $\gamma_{0}$ to each $\mu^{i}$, and they can be averaged to yield $\gamma_{1}$. If $|S|=M^{N}$ (that is, the collection $\{x^{i}_{m}\}$ satisfies a general-position-type condition), then $\gamma_{1}$ will be concentrated on $M$ points as well, and one may carry out further iterations. A conceptual problem with this application is that the Fréchet functional is not differentiable at discrete measures, so Algorithm 1 can no longer be viewed as gradient descent (but can still be seen as Procrustes averaging). Also, the Fréchet mean itself may fail to be unique. In simulations we observed very rapid convergence of this iteration to a Karcher mean, but the specific limit depended quite heavily on the initial point, and was usually not a Fréchet mean. For problems of moderate size, one can recast the problem of minimising the Fréchet functional as a linear program [9] and find an exact Fréchet mean. In fact, Anderes et al. [9] treat the more general problem where the measures are supported on a different number of points and are not constrained to be uniform on their supports.

An important issue more generally is that of efficient approximate numerical schemes for calculating Fréchet means in Wasserstein space. This is a very active field of research with a rapidly-growing literature (both in numerical analysis and in computer science), and a detailed survey is far beyond the scope of this paper. If one is content with an approximate solution, then there are several approaches suggested in the literature. Indicatively, let us mention Bonneel et al. [19] who use a tomographic perspective to reduce the problem to 1-dimensional computations; Carlier, Oberman & Oudet [22] who use nonsmooth optimisation techniques to solve a discretised version of the dual problem; Oberman & Ruan [59] exploit the sparsity of optimal plans to reduce the size of the linear program to a tractable one.

Another line of research involves entropic regularisation, where one adds an entropy term to the definition of the Wasserstein distance. This leads to a strictly convex problem that is far better behaved than the original problem. Though its solution no longer yields the actual mean, it can be thought of as a regularised surrogate Fréchet mean. In this direction, Cuturi & Doucet [27] employ differentiability properties and carry out what could be thought of as a “gradient descent”, a discrete analogue of Algorithm 1; Benamou et al.  [11] exploit the structure of the constraints as an intersection of convex sets by means of iterating Bregman projections that can be evaluated efficiently. Solomon et al. [70] extend this idea to the manifold setup, by convoluting with a heat kernel; and Cuturi & Peyré [28] employ the regularisation at the level of the dual, rather than the primal, problem. Recently, Rolet, Cuturi & Peyré [67] employed this technique in the context of dictionary learning; and Bonneel, Peyré & Cuturi [18] define a sort of “barycentric convex hull” of given histograms and show how to project a new histogram onto that convex hull.

Acknowledgements

This research was supported by a European Research Council Starting Grant Award to Victor M. Panaretos. Part of this work grew out of work presented at the Mathematical Biosciences Institute (Ohio State University), during the \hrefhttp://mbi.osu.edu/event/?id=162“Statistics of Time Warping and Phase Variation” Workshop, November 2012. We wish to acknowledge the stimulating environment offered by the Institute. We wish to warmly thank Prof. Clément Hongler for several useful discussions. We are also very thankful to two reviewers and an associate editor for their detailed and constructive feedback.

8 Supplementary Material

This section contains material supplementing the main article. The first section contains the proof that no further requirement except for finite second moments is needed for the convergence of the algorithm presented in the article. Next, we provide further details and theoretical results pertaining to the simulation scenarios described in Section 6. Finally, we provide all the proofs not included in the main body for tidiness, as well as additional technical details.

A complete proof of Lemma 4

In this section we show that condition (5.3) is not needed for (5.2) to hold. The idea is that (5.2) only requires a tiny bit more than finite second moments, and that is provided in Lemma 8.2. Throughout this section, all functions are assumed nonnegative (possibly infinite-valued) and defined on $[0,\infty)$ unless explicitly stated otherwise. We write $f(x)\in\omega(g(x))$ or $f\in\omega(g)$ if $f(x)/g(x)\to\infty$ as $x\to\infty$.

Lemma 8.1**.**

Let $f$ be integrable. Then there exists a continuous nondecreasing function $g\in\omega(1)$ such that $fg$ is integrable.

Proof.

Set $F(x)={\int_{x}^{\infty}\!f(t)\,\mathrm{d}t}$ and $g(x)=[F(x)]^{-1/2}$. Then a change of variables gives

[TABLE]

and $g(x)\to\infty$ because $F(x)\to 0$ as $x\to\infty$ by dominated convergence. ∎

Lemma 8.2**.**

Let $X$ be a random variable with $\mathbb{E}X^{2}<\infty$. Then there exists a convex nondecreasing function $H(x)\in\omega(x^{2})$ such that $\mathbb{E}H(X)<\infty$.

Proof.

Since

[TABLE]

there exists a function $g$ as in Lemma 8.1 such that

[TABLE]

where $G$ is the primitive of $g$ and $G(0)=0$. The properties of $g$ imply that $G$ is convex and invertible, and that for $y<x$,

[TABLE]

which, combined with $g(y)\to\infty$ as $y\to\infty$, yields

[TABLE]

so that $G(x)\in\omega(x)$. The function $H(x)=G(x^{2})$ then has all the desired properties. ∎

Proposition 8.3**.**

Equation (5.2) holds if merely

[TABLE]

Proof.

Let $X_{i}=\|Z^{i}\|$ where $Z^{i}\sim\mu^{i}$. Then there exist functions $g^{i}$ as in Lemma 8.1 with

[TABLE]

The same holds with $g^{i}$ replaced by $g=\min_{i}g^{i}$, which is still continuous, nondecreasing and divergent. Setting $H$ as in Lemma 8.2, we see that $H(x)\in\omega(x^{2})$ and

[TABLE]

Convexity of $H$ and $\|\cdot\|$ combined with monotonicity of $H$ yield

[TABLE]

where $M=\sum_{i=1}^{N}M^{i}/N$. This implies that for any $R>0$ and any $j>0$,

[TABLE]

and (5.2) follows because $H(y)\in\omega(y^{2})$. ∎

Details for the illustrative examples in Section 6

In this section we provide further details for finding the optimal maps in the examples of Section 6 and theoretical results about the Fréchet mean and the behaviour of the algorithm. Throughout this section, $\mu^{1},\dots,\mu^{N}$ are given measures and $\gamma_{0}$ is the initial point of Algorithm 1. We begin with two lemmas regarding compatibility of the measures as defined in Section 6.

Lemma 8.4** (Compatibility and Convergence).**

If $\mathbf{t}_{\mu^{1}}^{\mu^{i}}\circ\mathbf{t}_{\gamma_{0}}^{\mu^{1}}=\mathbf{t}_{\gamma_{0}}^{\mu^{i}}$ and $\mathbf{t}_{\mu^{1}}^{\mu^{j}}\circ\mathbf{t}_{\mu^{i}}^{\mu^{1}}=\mathbf{t}_{\mu^{i}}^{\mu^{j}}$ (in the relevant $L^{2}$ spaces) for all $i$ and all $j$, then Algorithm 1 converges after a single step.

Proof.

For all $i$, $j$ and $k$ we have $\mathbf{t}_{\mu^{j}}^{\mu^{k}}\circ\mathbf{t}_{\mu^{i}}^{\mu^{j}}=\mathbf{t}_{\mu^{i}}^{\mu^{k}}$, so that the optimal maps are admissible, and

[TABLE]

Boissard et al. [16] show that this is indeed the Fréchet mean. ∎

When $d=1$, all (diffuse) measures are compatible with each other, and Algorithm 1 converges after one step. Generally, the algorithm requires the calculation of $N$ pairwise optimal maps, and this can be reduced to $N-1$ if $\gamma_{0}=\mu^{1}$. This is the same computational complexity as the calculation of the iterated barycentre proposed in [16].

Measures on $\mathbb{R}^{d}$ that have a common dependence structure are compatible with each other. More precisely, we say that $C:[0,1]^{d}\to[0,1]$ is a copula if there exists a random vector $U$ with $U[0,1]$ margins and such that

[TABLE]

In other words, a copula is the restriction to $[0,1]^{d}$ of the probability distribution function of some $d$-dimensional random variable with uniform margins. See, for example, Nelsen [58] for an overview. Given a measure $\mu$ on $\mathbb{R}^{d}$ with distribution function $G$ and marginal distribution functions $G_{j}$, the copula associated with $\mu$ is a copula such that

[TABLE]

This equation defines $C$ uniquely if each marginal $G_{i}$ is continuous, which we shall assume for simplicity. (If some $G_{i}$ is discontinuous then $C$ might not be unique, but it always exists, see [58, Chapter 2].)

Lemma 8.5** (Compatibility and Copulae).**

Let $\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be regular. Then $\mu$ and $\nu$ have the same associated copula if and only if $\mathbf{t}_{\mu}^{\nu}$ takes the separable form

[TABLE]

The result can be obtained as a corollary of Cuesta-Albertos et al. [26, Theorem 2.9], but here is an alternative direct proof.

Proof.

If $\mu$ and $\nu$ have the same copula then

[TABLE]

where $G_{j}^{-1}(u_{j})$ is any number satisfying $G_{j}(G_{j}^{-1}(u_{j}))=u_{j}$ (such numbers exist because $G_{j}$ is surjective), and similarly for $F_{j}^{-1}$. Consequently, $F(x_{1},\dots,x_{d})=G(T_{1}(x_{1}),\dots,T_{d}(x_{d}))$ with $T_{j}=G_{j}^{-1}\circ F_{j}$. It follows that $\nu=(T_{1},\dots,T_{d})\#\mu$, and this map is optimal, hence equals $\mathbf{t}_{\mu}^{\nu}$, because the $T_{j}$’s are nondecreasing.

One proves the converse implication similarly: if $\mathbf{t}_{\mu}^{\nu}$ takes this form, then each $T_{j}$ needs to be nondecreasing. Since it must push $F_{j}$ forward to $G_{j}$, we have $T_{j}=G_{j}^{-1}\circ F_{j}$, and this yields the above equality for the copula. ∎

It is easy to see that if the optimal maps between each $\mu^{i}$ and each $\mu^{j}$ are of the form (8.1), then $\{\mu^{i}\}$ are compatible with other. This follows from this property holding for each marginal, and the possibility of working with the marginals separately; it has already been observed by Boissard et al. [16, Proposition 4.1]. This explains why the algorithm converges in one iteration for the example with the Frank copula.

Next, we give a convergence analysis for the Gaussian example.

Theorem 8.6** (Convergence in Gaussian case).**

Let $\mu^{i}\sim\mathcal{N}(0,S_{i})$ for $S_{i}$ positive definite, and let the initial point $\gamma_{0}=\mathcal{N}(0,\Gamma_{0})$ for positive definite $\Gamma_{0}$. Then the sequence of iterates generated by Algorithm 1 converges to the unique Fréchet mean of $(\mu^{1},\dots,\mu^{N})$.

Proof.

We first observe that for any centred measure $\mu$ with covariance matrix $S$,

[TABLE]

where $\delta_{0}$ is a dirac mass at the origin. (This follows from the singular value decomposition of $S$.) Next, each iteration stays (centred) Gaussian, say $\mathcal{N}(0,\Gamma_{k})$, because the optimal maps are linear; and since the iterates are absolutely continuous (Lemma 1), each $\Gamma_{k}$ is nonsingular.

Proposition 4 implies that $\mathrm{det}\Gamma_{k}$ is bounded below uniformly; on the other hand,

[TABLE]

is bounded uniformly, because $\{\gamma_{k}\}$ stays in a Wasserstein-compact set by Lemma 4. Let $C_{1}=\inf_{k}\mathrm{det}\Gamma_{k}>0$ and $C_{2}=\sup_{k}\mathrm{tr}\Gamma_{k}<\infty$. Then each eigenvalue $\lambda$ of $\Gamma_{k}$ is nonnegative, bounded above by $C_{2}$, and satisfies

[TABLE]

The matrices $\Gamma_{k}$ stay in a bounded set, and each limit point $\Gamma$ is positive definite because $x^{t}\Gamma x\geq C_{3}\|x\|^{2}$ for all $x\in\mathbb{R}^{d}$. Each limit point $\gamma$ of $\gamma_{k}$ is a Karcher mean by Theorem 3, and the limit must follow a $\mathcal{N}(0,\Gamma)$ distribution with $\Gamma$ (nonsingular) limit point of $\Gamma_{k}$ (e.g., by Lehmann–Scheffé’s theorem). Since $F^{\prime}(\gamma)=0$ everywhere on $\mathbb{R}^{d}$, $\gamma$ is the Fréchet mean by the discussion after Corollary 1. Every limit of $\gamma_{k}$ is the Fréchet mean and the sequence is compact, so $\gamma_{k}$ must converge to the Fréchet mean. ∎

Remark 1**.**

During the review process, a referee asked whether the result in Theorem 8.6 is related to the iteration $\Sigma\mapsto N^{-1}\sum_{i=1}^{N}(\Sigma^{1/2}\Sigma_{i}\Sigma^{1/2})^{1/2}$, introduced by Knott & Smith [48] and later considered in Xia et al. [76] for the Gaussian case. That iteration, however, is distinctly different from Algorithm 1, which Theorem 8.6 concerns (for instance, the latter involves inversion operations, which the former does not). As pointed out by Rüschendorf & Uckelmann [68, p. 6], the scheme $\Sigma\mapsto N^{-1}\sum_{i=1}^{N}(\Sigma^{1/2}\Sigma_{i}\Sigma^{1/2})^{1/2}$ is not known to converge, and indeed Theorem 8.6 does not furnish any additional insight on the matter.

In order to deal with the last example of Section 6, we need two more results. The first involves coupling measures of dimensions greater than one, while the second shows the equivariance of the Fréchet mean with respect to rotations.

Invoking the independence copula $C(u_{1},\dots,u_{d})=u_{1}\dots u_{d}$, a special case of Lemma 8.5 above is when the marginals of $\mu$ and $\nu$ are independent. In this independence case, it is possible in fact to replace the marginals by measures of arbitrary dimension:

Lemma 8.7**.**

Let $\mu^{1},\dots,\mu^{N}$ and $\nu^{1},\dots,\nu^{N}$ be regular measures in $\mathcal{P}_{2}(\mathbb{R}^{d_{1}})$ and $\mathcal{P}_{2}(\mathbb{R}^{d_{2}})$ with (unique) Fréchet means $\mu$ and $\nu$ respectively. Then the independent coupling $\mu\otimes\nu$ is the Fréchet mean of $\mu^{1}\otimes\nu^{1},\dots,\mu^{N}\otimes\nu^{N}$.

By induction (or a straightforward modification of the proof), one can show that the Fréchet mean of $(\mu^{i}\otimes\nu^{i}\otimes\rho^{i})$ is $\mu\otimes\nu\otimes\rho$, and so on. While we are confident this result should already be known, we could not find a reference, and thus we provide a full proof for completeness.

Proof.

Agueh & Carlier [2, Proposition 3.8] show that there exist convex lower semicontinuous potentials $\psi_{i}$ on $\mathbb{R}^{d_{1}}$ and $\varphi_{i}$ on $\mathbb{R}^{d_{2}}$ whose gradients push $\mu$ forward to $\mu^{i}$ and $\nu$ to $\nu^{i}$ respectively, and such that

[TABLE]

with equality $\mu$- and $\nu$-almost surely respectively. It is easy to see that the extensions $\tilde{\psi}_{i}(x,y)=\psi_{i}(x)$ and $\tilde{\varphi}_{i}(x,y)=\varphi_{i}(y)$ defined on $\mathbb{R}^{d_{1}+d_{2}}$ are convex lower semicontinuous functions whose sum $\phi_{i}$ is a convex function satisfying

[TABLE]

Clearly $\nabla\phi_{i}\#(\mu^{i}\otimes\nu^{i})=\mu\otimes\nu$ and

[TABLE]

with equality $\mu\otimes\nu$-almost surely. By the same Proposition 3.8 in [2], $\mu\otimes\nu$ is the Fréchet mean. ∎

Lemma 8.8**.**

If $\mu$ is the Fréchet mean of the regular measures $\mu^{1},\dots,\mu^{N}$, one with bounded density, and $U$ is orthogonal, then $U\#\mu$ is the Fréchet mean of $U\#\mu^{1},\dots,U\#\mu^{N}$.

Proof.

Bonneel et al. sketch a proof of this statement in [19, Proposition 1], and it also appears implicitly in Boissard et al. [16, Proposition 4.1]; we give an alternative argument here.

If $x\mapsto\varphi(x)$ is convex, then $x\mapsto\varphi(U^{-1}x)$ is convex with gradient $U\nabla\varphi(U^{-1}x)$ at (almost all) $x$ and conjugate $x\mapsto\varphi^{*}(U^{-1}x)$. If $\varphi_{i}$ are convex potentials with $\nabla\varphi_{i}\#\mu=\mu^{i}$, then $\nabla(\varphi_{i}\circ U^{-1})$ pushes $U\#\mu$ forward to $U\#\mu^{i}$ and by [2, Proposition 3.8]

[TABLE]

with equality for $\mu$-almost any $x$. A change of variables $y=Ux$ shows that the set of points $y$ such that $\sum(\varphi_{i}\circ U^{-1})^{*}(y)<N\|y\|^{2}/2$ is $(U\#\mu)$-negligible, completing the proof. ∎

We apply these results in the context of the simulated example in Section 6. If $Y=(y_{1},y_{2},y_{3})\sim\mu^{i}$, then the random vector $(x_{1},x_{2},x_{3})=X=U^{-1}Y$ has joint density

[TABLE]

so the probability law of $X$ is $\rho^{i}\otimes\nu^{i}$ with $\rho^{i}$ centred Gaussian with covariance matrix $\Sigma^{i}$ and $\nu^{i}$ having density $f^{i}$ on $\mathbb{R}$. By Lemma 8.7, the Fréchet mean of $(U^{-1}\#\mu^{i})$ is the product measure of that of $(\rho^{i})$ and that of $(\nu^{i})$; by Lemma 8.8, the Fréchet mean of $(\mu^{i})$ is therefore

[TABLE]

where $\Sigma$ is the Fréchet–Wasserstein mean of $\Sigma_{1},\dots,\Sigma_{N}$.

Starting at an initial point $\gamma_{0}=U\#(\mathcal{N}(0,\Sigma_{0})\otimes\nu_{0})$, with $\nu_{0}$ having continuous distribution $F_{\nu_{0}}$, the optimal maps are $U\circ\mathbf{t}_{0}^{i}\circ U^{-1}=\nabla(\varphi_{0}^{i}\circ U^{-1})$ with

[TABLE]

the gradients of the convex function

[TABLE]

where we identify $\mathbf{t}_{\gamma_{0}}^{\Sigma^{i}}$ with the positive definite matrix $(\Sigma^{i})^{1/2}[(\Sigma^{i})^{1/2}\Sigma_{0}(\Sigma^{i})^{1/2}]^{-1/2}(\Sigma^{i})^{1/2}$ that pushes forward $\mathcal{N}(0,\Sigma_{0})$ to $\mathcal{N}(0,\Sigma^{i})$. Due to the one-dimensionality, the algorithm finds the third component of the rotated measures after one step, but the convergence of the Gaussian component requires further iterations.

Proofs and Details Omitted from the Main Article

Proofs of statements from Section 3.1

Proof of Corollary 1, Section 3.2.3.

The characterisation of Karcher means follows immediately from Theorem 1. Now suppose that $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ is regular and $F^{\prime}(\mu)\neq 0\in L^{2}(\mu)$. The function $S=N^{-1}\sum_{i=1}^{N}\mathbf{t}_{\mu}^{\mu^{i}}$ is a gradient of a convex function and

[TABLE]

By assumption $W=S-\mathbf{i}\neq 0\in L^{2}(\mu)$. The measure $\nu_{s}=[\mathbf{i}+s(W-\mathbf{i})]\#\mu$ with $s\in(0,1)$ is such that $d(\nu_{s},\mu)=s\|W\|_{L^{2}(\mu)}$ and

[TABLE]

This means that when $s$ is small enough, $F(\nu_{s})<F(\mu)$, so $\mu$ cannot be the minimiser of $F$. Since $\bar{\mu}$ has to be regular [2, Proposition 5.1], necessity of $F^{\prime}(\bar{\mu})=0$ is proven. ∎

Proofs of statements from Section 4

Additional Details on the Proof of Theorem 5.

Write $M(\gamma)=\mathbb{E}[d^{2}(\Lambda,\gamma)]$. We wish to show that $M$ has a unique minimiser $\gamma$ and that $\gamma$ is supported on $K$. We first establish (weak) convexity of $M$. Indeed, for given measures $\gamma$ and $\rho$ and $0<t<1$,

[TABLE]

where $\pi_{\omega,\gamma}$ is the optimal coupling between $\Lambda=\Lambda(\omega)$ and $\gamma$. The measure $t\pi_{\omega,\gamma}+(1-t)\pi_{\omega,\rho}$ is a coupling between $\Lambda$ and $t\gamma+(1-t)\rho$, and this shows that $M_{\omega}$ is convex without any regularity assumptions on $\Lambda$. To upgrade to strict convexity when $\Lambda$ is regular, observe firstly that $M$ is finite on the set of probability measures supported on $K$. If $\Lambda$ is regular, then optimal measures are supported on graphs of functions:

[TABLE]

The measure $t\pi_{\omega,\gamma}+(1-t)\pi_{\omega,\rho}$ is supported on the graph of two functions, $T_{1}$ and $T_{2}$. It can only be optimal if it is supported on the graph of one function, and this will only happen if $T_{1}=T_{2}$, $\Lambda$-almost surely, that is, if $\gamma=\rho$. (See [5, Corollary 2.9] for a rigorous proof.) We can thus conclude that

[TABLE]

Since $M$ was already shown to be weakly convex in any case, it follows that

[TABLE]

Now we turn to the existence of a solution (once existence is established, uniqueness will follow from strict convexity). Let $\mathrm{proj}_{K}:\mathbb{R}^{d}\to K$ denote the projection onto the set $K$, which is well-defined since $K$ is closed and convex, and of course satisfies

[TABLE]

Since $\Lambda$ is concentrated on $K$, the above inequality holds $\Lambda$-almost surely with respect to $x$. Let $T$ be the optimal map from $\Lambda$ to $\gamma$ (a proper map almost surely, as argued above). Observe that

[TABLE]

since $(\mathrm{proj}_{K}\circ T)\#\Lambda=\mathrm{proj}_{K}\#(T\#\Lambda)=\mathrm{proj}_{K}\#\gamma$. This measure is concentrated on $K$, and taking expectations gives $M(\gamma)\geq M(\mathrm{proj}_{K}\#\gamma)$. Hence, the infimum of $M$ equals the infimum of $M$ on $\mathcal{P}(K)$, the collection of probability measures supported on $K$ (or else, we could project all the remaining mass to $K$ to reduce the total cost further). The restriction of $M$ to $\mathcal{P}(K)$ is a continuous functional on a compact set (measures whose support is contained in a common compactum are a compact set in Wasserstein space), and existence follows.

In order to establish that $\lambda$ minimises $M$, we need to justify the following facts:

[TABLE]

The space

[TABLE]

(endowed with the supremum norm $\|f\|_{\infty}=\sup_{x\in K}\|f(x)\|$) is a separable Banach space and therefore any random element $T:(\Omega,\mathcal{F},\mathbb{P})\to C_{b}(K,\mathbb{R}^{d})$ is Bochner measurable. Clearly

[TABLE]

so that $\Lambda=T\#\lambda$, when viewed as a random measure in the Wasserstein space, is a continuous function of $T$, and hence measurable. Again by continuity, $d^{2}(\Lambda,\theta):(\Omega,\mathcal{F},\mathbb{P})\to\mathbb{R}$ is measurable.

Let us now show the remaining two Fubini-type assertions. For simplicity, we assume that $K$ includes the origin. Then the convex potential $\phi$ of $T$ can be recovered as the line integral

[TABLE]

and the Legendre transform $\phi^{*}:K\to\mathbb{R}$ of $\phi$ by

[TABLE]

One can then verify the following properties (invoking the uniform continuity of $T$), which imply in particular that $\phi$ and $\phi^{*}$ are measurable:

$\phi$ is a continuous and bounded, hence an element of $C_{b}(K)$; 2. 2.

the same holds for $\phi^{*}$; 3. 3.

the map $T\mapsto\phi$ from $C_{b}(K,\mathbb{R}^{d})$ to $C_{b}(K)$ is Lipschitz; 4. 4.

the map $\phi\mapsto\phi^{*}$ from $C_{b}(K)$ to itself is Lipschitz.

Indeed, we write

[TABLE]

and notice that the first integral vanishes as $y\to x$ by uniform continuity. The last integral also vanishes by the Cauchy–Schwarz inequality, because $T$ is bounded (as a continuous function on $K$).

Again by the Cauchy–Schwarz inequality,

[TABLE]

so $\phi:C_{b}(K,\mathbb{R}^{d})\to C_{b}(K)$ is Lipschitz with constant $d_{K}=\sup_{x\in K}\|x\|$. It is also obvious that $\|\phi^{*}_{T}-\phi^{*}_{R}\|_{\infty}\leq\|\phi_{T}-\phi_{R}\|_{\infty}$ and that $\phi^{*}$ is bounded on $K$ because $\phi$ is bounded and $K$ are bounded. Uniform continuity can be verified directly as follows. Fix $\delta>0$ and $y,z\in K$ with $\|z-y\|<\delta$. Then for any $\epsilon>0$ we can pick some $x\in K$ such that

[TABLE]

Letting $\epsilon\to 0$ and since the supremum is finite, we see that $\phi^{*}(z)-\phi^{*}(y)\leq d_{K}\|y-z\|$; interchanging the roles of $y$ and $z$ above shows that in fact $|\phi^{*}(z)-\phi^{*}(y)|\leq d_{K}\|y-z\|$, so $\phi^{*}$ is even Lipschitz.

It now remains to show that

[TABLE]

Both equalities hold when $T$ is a simple function. Since $C_{b}(K,\mathbb{R}^{d})$ is separable, any $T$ can be approximated by simple functions $T_{n}$ such that $\|T_{n}\|\leq 2\|T\|$ almost surely. The assumption that $T$ takes values in $K$ almost surely implies that $\mathbb{E}T_{n}\to\mathbb{E}T$ in the Bochner sense, which means that $\|\mathbb{E}T_{n}-\mathbb{E}T\|_{\infty}\to 0$. Let $\phi_{n}$ be the convex potential of $T_{n}$. Then $\|\phi_{n}-\phi\|_{\infty}\leq\|T_{n}-T\|_{\infty}d_{K}$ and $\|\phi_{n}\|_{\infty}\leq\|T_{n}\|_{\infty}d_{K}$, which is integrable. It follows that $\mathbb{E}\phi_{n}\to\mathbb{E}\phi$ in the Bochner sense in $C_{b}(K)$, and in particular $\mathbb{E}\phi_{n}(x)\to\mathbb{E}\phi(x)$ for all $x\in K$, proving the first equality. The second equality is proven by a similar approximation argument. ∎

Proof of Lemma 5.

It is assumed that $\psi(z)=\psi_{1}(\|z\|)$ with $\psi_{1}$ non-increasing, strictly positive and

[TABLE]

Let $\Psi(A)={\int_{A}\!\psi(x)\,\mathrm{d}x}$ be the corresponding probability measure and recall that $\psi_{\sigma}(x)=\sigma^{-d}\psi(x/\sigma)$ for $\sigma>0$.

For $y\in K$ set $\tilde{\mu}_{y}=\delta\{y\}\ast\psi_{\sigma}$ and its restricted renormalized version $\mu_{y}=(1/\tilde{\mu}_{y}(K))\tilde{\mu}_{y}|_{K}$, so that $\widehat{\Lambda}_{i}=(1/m)\sum_{j=1}^{m}\mu_{x_{j}}$, and it is assumed that $m\geq 1$ and $x_{j}\in K$ (because $\Lambda_{i}(K)=1$).

One way (certainly not optimal, unless $m=1$) to couple $\widehat{\Lambda}_{i}$ with $\widetilde{\Pi}_{i}/m$ is to send the $1/m$ mass of $\mu_{x_{j}}$ to $x_{j}$. This gives

[TABLE]

However, for an arbitrary $y\in K$,

[TABLE]

The last displayed integral is bounded by 1. Hence, we seek a lower bound, uniformly in $y$ and $\sigma$, for

[TABLE]

Since $K-y$ is a convex set that contains the origin, the collection of sets $\{\sigma^{-1}(K-y)\}_{\sigma>0}$ is increasing as $\sigma\searrow 0$. Consequently, if $\sigma\leq 1$,

[TABLE]

Here $d_{K}=\sup\{\|x-y\|:x,y\in K\}$ is the (finite) diameter of $K$, and we have used the monotonicity of $\psi_{1}$.

It follows that for $C_{\psi,K}=[\psi_{1}(d_{K})\mathrm{Leb}(K)]^{-1}<\infty$ (depending only on $\psi$ and $K$),

[TABLE]

Since the bound is uniform in $y$, the proof is complete. In the context of Remark 5, one simply needs to replace the term $\psi_{1}(d_{K})$ in $C_{\psi,K}$ by $\delta_{\psi}(d_{K})$. ∎

Proof of Theorem 6.

Convergence in probability in part (1) follows as in [61, pp. 793–794], using (5.6). For convergence almost surely, let $a=(a_{1},\dots,a_{d})\in\mathbb{R}^{d}$. A straightforward generalisation of the argument in [61, pp. 794–795] gives

[TABLE]

where for $-\infty\leq a_{i}\leq b_{i}\leq\infty$ ($i=1,\dots,d$) we denote

[TABLE]

Consequently

[TABLE]

For a general $a\in\mathbb{R}^{d}$ there exist sequences $a^{k}\nearrow a\swarrow b^{k}$ with $a^{k},b^{k}\in\mathbb{Q}^{d}$ (that is, $a_{i}^{k}\nearrow a_{i}\swarrow b_{i}^{k}$ for any coordinate $i$). Since for any $k$

[TABLE]

it follows that with probability one

[TABLE]

as the sequence of sets at the right-hand side converges monotonically to the empty set.

Similarly, with probability one

[TABLE]

and the right-hand side vanishes because $\Lambda_{i}$ is assumed absolutely continuous (the set $(-\infty,a]\setminus(-\infty,a)$ is union of $d$ $d-1$-dimensional rays). Specifying $a=\infty$ shows that almost surely $\widetilde{\Pi}_{i}(K)/\tau_{n}\to 1$ and we conclude that almost surely $\widetilde{\Pi}_{i}/\widetilde{\Pi}_{i}(K)\to\Lambda_{i}$ weakly. Further, $d(\widehat{\Lambda}_{i},\widetilde{\Pi}_{i}/N_{i})\to 0$ since $\sigma_{n}\to 0$ by Lemma 5.

We sketch the main ideas of the proof of (2); more details can be found in [61, pp. 795–797]. We wish to show that

[TABLE]

In order to do this we write

[TABLE]

where we introduce the empirical Fréchet functional

[TABLE]

Since for any three probability measures on $K$ it holds that

[TABLE]

we see that

[TABLE]

Each $X_{ni}$ is a function of $T_{i}$, $\Pi_{i}^{(n)}$ and $\sigma_{i}^{(n)}$, and $0\leq X_{ni}\leq d_{K}$. If $\sigma_{i}^{(n)}$ is a function of $\widetilde{\Pi}_{i}^{(n)}=T_{i}\#\Pi_{i}^{(n)}$ only, then $X_{ni}$ are iid across $i$. Part (1) shows that $X_{n1}\to 0$ in probability and by the bounded convergence theorem $\mathbb{E}\overline{X}_{n}=\mathbb{E}X_{n1}\stackrel{{\scriptstyle p}}{{\to}}0$ and therefore the above expression converges to 0 in probability. In general, $L^{1}$-convergence of random variables does not guarantee convergence almost surely. As we deal with averages, however, almost sure convergence can be established: let $Y_{ni}=X_{ni}-\mathbb{E}X_{ni}\in[-d_{K},d_{K}]$. Then $Y_{ni}$ are mean zero iid random variables, so that

[TABLE]

By the Borel–Cantelli lemma, $|\overline{X}_{n}-\mathbb{E}\overline{X}_{n}|\stackrel{{\scriptstyle as}}{{\to}}0$, hence $\overline{X}_{n}\stackrel{{\scriptstyle as}}{{\to}}0$.

If the smoothing is not carried out independently across trains, then $X_{ni}$ may be correlated across $i$. In that case, one can introduce the functional $M_{n}^{*}(\gamma)=n^{-1}\sum_{i=1}^{n}d^{2}\left(\widetilde{\Pi}_{i}/N_{i},\gamma\right)$ and proceed as in [61]. For $M_{n}^{*}$ to be well-defined one may use Lemma 6 and that requires $\tau_{n}/\log n\to\infty$.

Finally, observe that by the strong law of large numbers $M_{n}(\gamma)\stackrel{{\scriptstyle as}}{{\to}}M(\gamma)$ for all $\gamma\in P(K)$. That the convergence is uniform follows from the equicontinuity of the collection $\{M_{n}\}_{n=1}^{\infty}$ (they are $2d_{K}$-Lipschitz). We have thus established

[TABLE]

By standard arguments, the minimiser $\widehat{\lambda}_{n}$ of $\widehat{M}_{n}$ converges to the minimiser $\lambda^{*}$ of $M$, since the latter is unique by Theorem 5. But $\lambda^{*}=\lambda$ by the hypothesis. ∎

Proofs of statements from Section 5.5

Proof of Lemma 7.

For any $1>\epsilon>0$ there exists $0<t_{\epsilon}$ such that for $t<t_{\epsilon}$,

[TABLE]

Fix $z$ such that $t=t(z)=\|z-x_{0}\|<t_{\epsilon}$. The intersection of $B_{t}(x_{0})$ with $B_{2\epsilon t}(z)$ includes a ball of radius $\epsilon t$ centred at $y=x_{0}+(1-\epsilon)(z-x_{0})$, so that

[TABLE]

It follows that $G\cap B_{2\epsilon t}(z)$ is nonempty. In other words: for any $\epsilon>0$ there exists $t_{\epsilon}$ such that if $\|z-x_{0}\|<t_{\epsilon}$, then there exists $x\in G$ with $\|z-x\|\leq 2\epsilon t(z)=2\epsilon\|z-x_{0}\|$. This means precisely that $\delta(z)=o(\|z-x_{0}\|)$ as $z\to x_{0}$. ∎

Proof of Lemma 12.

Assume $\epsilon<\rho$ (there is nothing to prove otherwise). Take a corner of the $\ell_{\infty}$ ball of radius $\rho^{\prime}<\rho$ around $x_{0}$,

[TABLE]

and write $y=\sum a_{i}z_{i}$ as a (finite) convex combination of elements of $Z$. Then $\tilde{y}=\sum a_{i}\tilde{z}_{i}\in\mathrm{conv}(\tilde{Z})$ is such that $\|\tilde{y}-y\|_{\infty}\leq\epsilon$. It follows that $\tilde{y}$ lies at the same quadrant as $y$ with each coordinate larger in absolute value than $\rho^{\prime}-\epsilon$. In other words, $\tilde{y}$ is “more extreme” than the corner

[TABLE]

of the $\ell_{\infty}$-ball $B_{\rho^{\prime}-\epsilon}^{\infty}(x_{0})$. Since this is true for all the corners, $\mathrm{conv}(\tilde{Z})\supseteq B_{\rho^{\prime}-\epsilon}(x_{0})$ for any $\rho^{\prime}<\rho$. Now let $\rho^{\prime}\nearrow\rho$ to conclude. ∎