© The Author(s) 2020

V. M. Panaretos, Y. ZemelAn Invitation to Statistics in Wasserstein SpaceSpringerBriefs in Probability and Mathematical Statisticshttps://doi.org/10.1007/978-3-030-38438-8_4

4. Phase Variation and Fréchet Means

Victor M. Panaretos¹  and Yoav Zemel²

(1)

Institute of Mathematics, EPFL, Lausanne, Switzerland

(2)

Statistical Laboratory, University of Cambridge, Cambridge, UK

 

Electronic Supplementary Material The online version of this chapter (https://​doi.​org/​10.​1007/​978-3-030-38438-8_​4) contains supplementary material.

For instance, the arithmetic average of two scalar Gaussians N(μ ₁, 1) and N(μ ₂, 1) will be their mixture with equal weights, but their Fréchet–Wasserstein average will be the Gaussian (see Lemma 4.2.1), which is arguably more representative from an intuitive point of view. In much the same way, the Fréchet–Wasserstein average of probability measures representing some type of object (e.g., normalised greyscale images of faces) will also be an object of the same type. This sort of phenomenon is well-known in manifold statistics, more generally, and is arguably one of the key motivations to account for the non-linear geometry of the sample space, rather than imbed it into a larger linear space and use the addition operation.

Why is it relevant to construct the Fréchet mean of a collection of measures with respect to the Wasserstein metric? A simple answer is that this kind of average will often express a more natural notion of “typical” realisation of a random probability distribution than an arithmetic average.¹ Much more can be said, however, in that the Wasserstein–Fréchet mean and the closely related notion of an optimal multicoupling arise canonically as the appropriate framework for the formulation and solution to the problem of separation of amplitude and phase variation of a point process. It would almost seem that Wasserstein–Fréchet means were “made” for precisely this problem.

When analysing the (co)variation of a real-valued stochastic process {Y (x) : x ∈ K} over a convex compact domain K, it can be broadly said that one may distinguish two layers of variation:

• Amplitude variation. This is the “classical” variation that one would also encounter in multivariate analysis, and refers to the stochastic fluctuations around a mean level, usually encoded in the covariance kernel, at least up to second order.

  In short, this is variation “in the y-axis” (ordinate).

• Phase variation. This is a second layer of non-linear variation peculiar to continuous domain stochastic processes, and is rarely—if ever—encountered in multivariate analysis. It arises as the result of random changes (or deformations) in the time scale (or the spatial domain) of definition of the process. It can be conceptualised as a composition of the stochastic process with a random transformation (warp map) acting on its domain.

  This is variation “in the x-axis” (abscissa).

The terminology on amplitude/phase variation is adapted from random trigonometric functions, which may vary in amplitude (oscillations in the range of the function) or phase (oscillations in the domain of the function). Failing to properly account for the superposition of these two forms of variation may entirely distort the findings of a statistical analysis of the random function (see Sect. 4.1.1). Consequently, it is an important problem to be able to separate the two, thus correctly accounting for the distinct contribution of each. The problem of separation is also known as that of registration (Ramsay and Li [108]), synchronisation (Wang and Gasser [129]), or multireference alignment (Bandeira et al. [16]), though in some cases these terms refer to a simpler problem where there is no amplitude variation at all.

Phase variation naturally arises in the study of random phenomena where there is no absolute notion of time or space, but every realisation of the phenomenon evolves according to a time scale that is intrinsic to the phenomenon itself, and (unfortunately) unobservable. Processes related to physiological measurements, such as growth curves and neuronal signals, are usual suspects. Growth curves can be modelled as continuous random functions (functional data), whereas neuronal signals are better modelled as discrete random measures (point processes). We first describe amplitude/phase variation in the former² case, as that is easier to appreciate, before moving on to the latter case, which is the main subject of this chapter.

4.1 Amplitude and Phase Variation

4.1.1 The Functional Case

Let K denote the unit cube . A real random function Y = (Y (x) : x ∈ K) can, broadly speaking, have two types of variation. The first, amplitude variation, results from Y (x) being a random variable for every x and describes its fluctuations around the mean level , usually encoded by the variance varY (x). For this reason, it can be referred to as “variation in the y-axis”. More generally, for any finite set x ₁, …, x _n, the n × n covariance matrix with entries κ(x _i, x _j) = cov[Y (x _i), Y (x _j)] encapsulates (up to second order) the stochastic deviations of the random vector (Y (x ₁), …, Y (x _n)) from its mean, in analogy with the multivariate case. Heuristically, one then views amplitude variation as the collection κ(x, y) for x, y ∈ K in a sense we discuss next.

One typically views Y as a random element in the separable Hilbert space L ₂(K), assumed to have and continuous sample paths, so that in particular Y (x) is a random variable for all x ∈ K. Then the mean function

and the covariance kernel

are well-defined and finite; we shall assume that they are continuous, which is equivalent to Y being mean-square continuous:

The covariance kernel κ gives rise to the covariance operator , defined by

a self-adjoint positive semidefinite Hilbert–Schmidt operator on L ₂(K). The justification to this terminology is the observation that when m = 0, for all bounded f, g ∈ L ₂(K),

and so, without the restriction to m = 0,

The covariance operator admits an eigendecomposition such that , and (ϕ _k) is an orthonormal basis of L ₂(K). One then has the celebrated Karhunen–Loève expansion

A major feature in this expansion is the separation of the functional part from the stochastic part: the functions ϕ _k(x) are deterministic; the random variables ξ _k are scalars. This separation actually holds for any orthonormal basis; the role of choosing the eigenbasis of is making ξ _k uncorrelated:

vanishes when k≠l and equals r _k otherwise. For this reason, it is not surprising that using as ϕ _k the eigenfunctions yields the optimal representation of Y . Here, optimality is with respect to truncations: for any other basis (ψ _k) and any M,

so that (ϕ _k) provides the best finite-dimensional approximation to Y . The approximation error on the right-hand side equals

and depends on how quickly the eigenvalues of decay.

One carries out inference for m and κ on the basis of a sample Y ₁, …, Y _n by

and

from which one proceeds to estimate and its eigendecomposition.

We have seen that amplitude variation in the sense described above is linear and dealt with using linear operations. There is another, qualitatively different type of variation, phase variation, that is non-linear and does not have an obvious finite-dimensional analogue. It arises when in addition to the randomness in the values Y (x) itself, an extra layer of stochasticity is present in its domain of definition. In mathematical terms, there is a random invertible warp function (sometimes called deformation or warping) T : K → K and instead of Y (x), one observes realisations from

For this reason, phase variation can be viewed as “variation in the x-axis”. When d = 1, the set K is usually interpreted as a time interval, and then the model stipulates that each individual has its own time scale. Typically, the warp function is assumed to be a homeomorphism of K independent of Y and often some additional smoothness is imposed, say T ∈ C ². One of the classical examples is growth curves of children, of which a dataset from the Berkeley growth study (Jones and Bayley [73]) is shown in Fig. 4.1. The curves are the derivatives of the height of a sample of ten girls as a function of time, from birth until age 18. One clearly notices the presence of the two types of variation in the figure. The initial velocity for all children is the highest immediately or shortly after birth, and in most cases decreases sharply during the first 2 years. Then follows a period of acceleration for another year or so, and so on. Despite presenting qualitatively similar behaviour, the curves differ substantially not only in the magnitude of the peaks but also in their location. For instance, one red curve has a local minimum at the age of three, while a green one has a local maximum at almost that same time point. It is apparent that if one tries to estimate the mean function by averaging the curves at each time x, the shape of the resulting estimate would look very different from each of the curves. Thus, this pointwise averaging (known as the cross-sectional mean) fails to represent the typical behaviour. This phenomenon is seen more explicitly in the next example.

Fig. 4.1

Derivatives of growth curves of ten girls from the Berkeley dataset. The data and the code for the figure are from the R package fda (Ramsay et al. [111])

The terminology of amplitude and phase comes from trigonometric functions, from which we derive an artificial example that illustrates the difficulties of estimation in the presence of phase variation. Let A and B be symmetric random variables and consider the random function

(4.1)

(Strictly speaking, x↦x + B is not from [0, 1] to itself; for illustration purposes, we assume in this example that .) The random variable A generates the amplitude variation, while B represents the phase variation. In Fig. 4.2, we plot four realisations and the resulting empirical means for the two extreme scenarios where B = 0 (no phase variation) or A = 1 (no amplitude variation). In the left panel of the figure, we see that the sample mean (in thick blue) lies between the observations and has a similar form, so can be viewed as the curve representing the typical realisation of the random curve. This is in contrast to the right panel, where the mean is qualitatively different from all curves in the sample: though periodicity is still present, the peaks and troughs have been flattened, and the sample mean is much more diffuse than any of the observations.

Fig. 4.2

Four realisations of (4.1) with means in thick blue. Left: amplitude variation (B = 0); right: phase variation (A = 1)

The phenomenon illustrated in Fig. 4.2 is hardly surprising, since as mentioned earlier amplitude variation is linear while phase variation is not, and taking sample means is a linear operation. Let us see in formulae how this phenomenon occurs. When A = 1 we have

Since B is symmetric the second term vanishes, and unless B is trivial the expectation of the cosine is smaller than one in absolute value. Consequently, the expectation of is the original function multiplied by a constant of magnitude strictly less than one, resulting in peaks of smaller magnitude.

In the general case, where and Y and T are independent, we have

and

From this, several conclusions can be drawn. Let be the conditional mean function given T. Then the value of the mean function itself, , at x ₀ is determined not by a single point, say x, but rather by all the values of m at the possible outcomes of T ⁻¹(x). In particular, if x ₀ was a local maximum for m, then will typically be strictly smaller than m(x ₀); the phase variation results in smearing m.

At this point an important remark should be made. Whether or not phase variation is problematic depends on the specific application. If one is interested indeed in the mean and covariance functions of , then the standard empirical estimators will be consistent, since itself is a random function. But if it is rather m, the mean of Y , that is of interest, then the confounding of the amplitude and phase variation will lead to inconsistency. This can also be seen from the formula

The above series is not the Karhunen–Loève expansion of ; the simplest way to notice this is the observation that ϕ _k(T ⁻¹(x)) includes both the functional component ϕ _k and the random component T ⁻¹(x). The true Karhunen–Loève expansion of will in general be qualitatively very different from that of Y , not only in terms of the mean function but also in terms of the covariance operator and, consequently, its eigenfunctions and eigenvalues. As illustrated in the trigonometric example, the typical situation is that the mean is more diffuse than m, and the decay of the eigenvalues of the covariance operator is slower than that of r _k; as a result, one needs to truncate the sum at high threshold in order to capture a substantial enough part of the variability. In the toy example (4.1), the Karhunen–Loève expansion has a single term besides the mean if B = 0, while having two terms if A = 1.

When one is indeed interested in the mean m and the covariance κ, the random function T pertaining to the phase variation is a nuisance parameter. Given a sample , i = 1, …, n, there is no point in taking pointwise means of , because the curves are misaligned; should not be compared with , but rather with . To overcome this difficulty, one seeks estimators such that

is approximately Y _i(x). In other words, one tries to align the curves in the sample to have a common time scale. Such a procedure is called curve registration. Once registration has been carried out, one proceeds the analysis on assuming only amplitude variation is now present: estimate the mean m by

and the covariance κ by its analogous counterpart. Put differently, registering the curves amounts to separating the two types of variation. This step is crucial regardless of whether the warp functions are considered as nuisance or an analysis of the warp functions is of interest in the particular application.

There is an obvious identifiability problem in the model . If S is any (deterministic) invertible function, then the model with (Y, T) is statistically indistinguishable from the model with (Y ∘ S, T ∘ S). It is therefore often assumed that is the identity and in addition, in nearly all application, that T is monotonically increasing (if d = 1).

Discretely observed data. One cannot measure the height of person at every single instant of her life. In other words, it is rare in practice that one has access to the entire curve. A far more common situation is that one observes the curves discretely, i.e., at a finite number of points. The conceptually simplest setting is that one has access to a grid x ₁, …, x _J ∈ K, and the data come in the form

possibly with measurement error. The problem is to find, given , consistent estimators of T _i and of the original, aligned functions Y _i.

In the bibliographic notes, we review some methods for carrying out this separation of amplitude and phase variation. It is fair to say that no single registration method arises as the canonical solution to the functional registration problem. Indeed, most need to make additional structural and/or smoothness assumptions on the warp maps, further to the basic identifiability conditions requiring that T be increasing and that equal the identity. We will eventually see that, in contrast, the case of point processes (viewed as discretely observed random measures) admits a canonical framework, without needing additional assumptions.

4.1.2 The Point Process Case

A point process is the mathematical object that represents the intuitive notion of a random collection of points in a space . It is formally defined as a measurable map Π from a generic probability space into the space of (possibly infinite) Borel integer-valued measures of in such a way that Π(B) is a measurable real-valued random variable for all Borel subsets B of . The quantity Π(B) represents the random number of points observed in the set B. Among the plethora of books on point processes, let us mention Daley and Vere-Jones [41] and Karr [79]. Kallenberg [75] treats more general objects, random measures, of which point processes are a peculiar special case. We will assume for convenience that Π is a measure on a compact subset .

Amplitude variation of Π can be understood in analogy with the functional case. One defines the mean measure

and, provided that , the covariance measure

the latter being a finite signed Borel measure on K. Just like in the functional case, these two objects encapsulate the (second-order) amplitude variation³ properties of the law of Π.

Given a sample Π ₁, …, Π _n of independent point processes distributed as Π, the natural estimators

are consistent and the former asymptotically normal [79, Proposition 4.8].

Phase variation then pertains to a random warp function T : K → K (independent of Π) that deforms Π: if we denote the points of Π by x ₁, …, x _K (with K random), then instead of (x _i), one observes T(x ₁), …, T(x _K). In symbols, this means that the data arise as . We refer to Π as the original point processes, and as the warped point processes. An example of 30 warped and unwarped point processes is shown in Fig. 4.3. The point patterns in both panels present a qualitatively similar structure: there are two peaks of high concentration of points, while few points appear between these peaks. The difference between the two panels is in the position and concentration of those peaks. In the left panel, only amplitude variation is present, and the location/concentration of the peaks is the same across all observations. In contrast, phase variation results in shifting the peaks to different places for each of the observations, while also smearing or sharpening them. Clearly, estimation of the mean measure of a subset A by averaging the number of observed points in A would not be satisfactory as an estimator of λ when carried out with the warped data. As in the functional case, it will only be consistent for the measure defined by

and misses most (or at least a significant part) of the bimodal structure of λ and is far more diffuse.

Fig. 4.3

Unwarped (left) and warped Poisson point processes

Since Π and T are independent, the conditional expectation of given T is

Consequently, we refer to Λ = T#λ as the conditional mean measure. The problem of separation of amplitude and phase variation can now be stated as follows. On the basis of a sample , find estimators of (T _i) and (Π _i). Registering the point processes amounts to constructing estimators, registration maps , such that the aligned points

are close to the original points Π _i.

Remark 4.1.1 (Poisson Processes)

A special but important case is that of a Poisson process. Gaussian processes probably yield the most elegant and rich theory in functional data analysis, and so do Poisson processes when it comes to point processes. We say that Π is a Poisson process when the following two conditions hold. (1) For any disjoint collection (A ₁, …, A _n) of Borel sets, the random variables Π(A ₁), …, Π(A _n) are independent; and (2) for every Borel A ⊆ K, Π(A) follows a Poisson distribution with mean λ(A):

Conditional upon T, the random variables , k = 1, …, n are independent as the sets (T ⁻¹(A _k)) are disjoint, and follows a Poisson distribution with mean λ(T ⁻¹(A)) = Λ(A). This is precisely the definition of a Cox process : conditional upon the driving measure Λ, is a Poisson process with mean measure λ. For this reason, it is also called a doubly stochastic process ; in our context, the phase variation is associated with the stochasticity of Λ while the amplitude one is associated with the Poisson variation conditional upon Λ.

As in the functional case there are problems with identifiability: the model (Π, T) cannot be distinguished from the model (S#Π, T ∘ S ⁻¹) for any invertible S : K → K. It is thus natural to assume that is the identity map⁴ (otherwise set , i.e., replace Π by and T by ).

Constraining T to have mean identity is nevertheless not sufficient for the model to be identifiable. The reason is that given the two point sets and Π, there are many functions that push forward the latter to the former. This ambiguity can be dealt with by assuming some sort of regularity or parsimony for T. For example, when K = [a, b] is a subset of the real line, imposing T to be monotonically increasing guarantees its uniqueness. In multiple dimensions, there is no obvious analogue for increasing functions. One possible definition is the monotonicity described in Sect. 1.​7.​2:

This property is rather weak in a sense we describe now. Let and write y ≥ x if and only if y _i ≥ x _i for i = 1, 2. It is natural to expect the deformations to maintain the lexicographic order in :

If we require in addition that the ordering must be preserved for all quadrants: for z = T(x) and w = T(y)

then monotonicity is automatically satisfied. In that sense, it is arguably not very restrictive.

Monotonicity is weaker than cyclical monotonicity (see (1.​10) with y _i = T(x _i)), which is itself equivalent to the property of being the subgradient of a convex function. But if extra smoothness is present and T is a gradient of some function , then ϕ must be convex and T is then cyclically monotone. Consequently, we will make the following assumptions:

• the expected value of T is the identity;

• T is a gradient of a convex function.

In the functional case, at least on the real line, these two conditions are imposed on the warp functions in virtually all applications, often accompanied with additional assumptions about smoothness of T, its structural properties, or its distance from the identity. In the next section, we show how these two conditions alone lead to the Wasserstein geometry and open the door to consistent, fully nonparametric separation of the amplitude and phase variation.

4.2 Wasserstein Geometry and Phase Variation

4.2.1 Equivariance Properties of the Wasserstein Distance

A first hint to the relevance of Wasserstein metrics in for deformations of the space is that for all p ≥ 1 and all ,

where δ _x is as usual the Dirac measure at . This is in contrast to metrics such as the bounded Lipschitz distance (that metrises weak convergence) or the total variation distance on . Recall that these are defined by

so that

In words, the total variation metric “does not see the geometry” of the space . This is less so for the bounded Lipschitz distance that does take small distances into account but not large ones.

Another property (shared by BL and TV) is equivariance with respect to translations. It is more convenient to state it using the probabilistic formalism of Sect. 1.​2. Let X ∼ μ and Y ∼ ν be random elements in , a be a fixed point in , X′ = X + a and Y ′ = Y + a. Joint couplings Z′ = (X′, Y ′) are precisely those that take the form (a, a) + Z for a joint coupling Z = (X, Y ). Thus

where δ _a is a Dirac measure at a and ∗ denotes convolution.

This carries over to Fréchet means in an obvious way.

Lemma 4.2.1 (Fréchet Means and Translations)

Let Λ be a random measure in with finite Fréchet functional and . Then γ is a Fréchet mean of Λ if and only if γ ∗ δ _a is a Fréchet mean of Λ ∗ δ _a.

The result holds for other values of p, in the formulation sketched in the bibliographic notes of Chap. 2. In the quadratic case, one has a simple extension to the case where only one measure is translated. Denote the first moment (mean) of by

(When is infinite-dimensional, this can be defined as the unique element satisfying

By an equivalence of couplings similar to above, we obtain

which is minimised at a = m(μ) − m(ν). This leads to the following conclusion:

Proposition 4.2.2 (First Moment of Fréchet Mean)

Let Λ be a random measure in with finite Fréchet functional and Fréchet mean γ. Then

4.2.2 Canonicity of Wasserstein Distance in Measuring Phase Variation

The purpose of this subsection is to show that the standard functional data analysis assumptions on the warp function T, having mean identity and being increasing, are equivalent to purely geometric conditions on T and the conditional mean measure Λ = T#λ. Put differently, if one is willing to assume that and that T is increasing, then one is led unequivocally to the problem of estimation of Fréchet means in the Wasserstein space . When , “increasing” is interpreted as being the gradient of a convex function, as explained at the end of Sect. 4.1.2.

The total mass is invariant under the push-forward operation, and when it is finite, we may assume without loss of generality that it is equal to one, because all the relevant quantities scale with the total mass. Indeed, if λ = τμ with μ probability measure and τ > 0, then T#λ = τ × T#μ, and the Wasserstein distance (defined as the infimum-over-coupling integrated cost) between τμ and τν is τW _p(μ, ν) for μ, ν probabilities.

We begin with the one-dimensional case, where the explicit formulae allow for a more transparent argument, and for simplicity we will assume some regularity.

Assumptions 2

The domain is a nonempty compact convex set (an interval), and the continuous and injective random map (a random element in C _b(K)) satisfies the following two conditions:

1. (A1)

   Unbiasedness: for all x ∈ K.

    
2. (A2)

   Regularity: T is monotone increasing.

    

The relevance of the Wasserstein geometry to phase variation becomes clear in the following proposition that shows that Assumptions 2 are equivalent to geometric assumptions on the Wasserstein space .

Proposition 4.2.3 (Mean Identity Warp Functions and Fréchet Means in )

Let compact and convex and continuous. Then Assumptions 2hold if and only if, for any supported on K such that , the following two conditions are satisfied:

1. (B1)
   Unbiasedness: for any

   
    
2. (B2)
   Regularity: if is such that T#λ = Q#λ, then with probability one

   
    

These assumptions have a clear interpretation: (B1) stipulates that λ is a Fréchet mean of the random measure Λ = T#λ, while (B2) states that T must be the optimal map from λ to Λ, that is, .

Proof

If T satisfies (B2) then, as an optimal map, it must be nondecreasing λ-almost surely. Since λ is arbitrary, T must be nondecreasing on the entire domain K. Conversely, if T is nondecreasing, then it is optimal for any λ. Hence (A2) and (B2) are equivalent.

Assuming (A2), we now show that (A1) and (B1) are equivalent. Condition (B1) is equivalent to the assertion that for all ,

which is in turn equivalent to (see Sect. 3.​1.​4). Condition (A2) and the assumptions on T imply that F _Λ(x) = F _λ(T ⁻¹(x)). Suppose that F _λ is invertible (i.e., continuous and strictly increasing on K). Then . Thus (B1) is equivalent to for all x in the range of , which is K. The assertion that (A1) implies (B1), even if F _λ is not invertible, is proven in the next theorem (Theorem 4.2.4) in a more general context.

The situation in more than one dimension is similar but the proof is less transparent. To avoid compactness assumptions, we introduce the following power growth condition (taken from Agueh and Carlier [2]) of continuous functions that grow like ∥⋅∥^q (q ≥ 0):

with the norm . The space is defined similarly, with f taking values in instead of , and the norm will be denoted in the same way. These are nonseparable Banach spaces.

Theorem 4.2.4 (Mean Identity Warp Functions and Fréchet Means)

Fix and let be a (Bochner measurable) random optimal map with (Bochner) mean identity and such that . Then Λ = t#λ has Fréchet mean λ:

The generalisation with respect to the one-dimensional result is threefold. Firstly, since our main interest is the implication (A1–A2) ⇒ (B1–B2), we need not assume T to be injective. Secondly, the support of λ is not required to be compact. Lastly, the result holds in arbitrary dimension, including infinite-dimensional separable Hilbert spaces . In particular, if t is a linear map, then coincides with the operator norm of t, so the assumption is that t be a bounded self-adjoint nonnegative operator with mean identity and finite expected operator norm.

Proof

Optimality of t ensures that it has a convex potential ϕ, and strong and weak duality give

Formally taking expectations, using Fubini’s theorem and that (since is the identity) yields

as required. The rigorous mathematical justification for this is given on page 88 in the supplement.

Remark 4.2.5

The “natural” space fortwould be , but without the continuity assumption, the result may fail (Álvarez-Esteban et al. [ 9, Example 3.1]). A simple argument shows that the growth condition imposed by the G ₁assumption is minimal; see page 89 in the supplement or Galasso et al. [ 58].

Remark 4.2.6

The same statement holds if is replaced by a (Borel) convex subset K thereof. The integrals will then be taken on K, showing that λ minimises the Fréchet functional among measures supported on K, and, by continuity, on . By Proposition 3.2.4 , λ is a Fréchet mean.

4.3 Estimation of Fréchet Means

4.3.1 Oracle Case

In view of the canonicity of the Wasserstein geometry in Sect. 4.2.2, separation of amplitude and phase variation of the point processes essentially requires computing Fréchet means in the 2-Wasserstein space. It is both conceptually important and technically convenient to introduce the case where an oracle reveals the conditional mean measures Λ = T#λ entirely. Thus, assuming that is the unique Fréchet mean of a random measure Λ, the goal is to estimate the structural mean λ on the basis of independent and identically distributed realisations Λ ₁, …, Λ _n of λ.

Given that λ is defined as the minimiser of the Fréchet functional

it is natural to estimate λ by a minimiser, say λ _n, of the empirical Fréchet functional

A minimiser λ _n exists by Corollary 3.​1.​3. When , λ _n can be seen to be an unbiased estimator of λ in a generalised sense of Lehmann [88] (see Sect. 4.3.5).

The warp maps (and their inverses) can then be estimated as the optimal maps from λ _n to each Λ _i (see Sect. 4.3.4).

4.3.2 Discretely Observed Measures

In practice, one does not have the fortune of fully observing the inherently infinite-dimensional objects Λ ₁, …, Λ _n. A far more realistic scenario is that one only has access to a discrete version of Λ _i, say . The simplest situation is when arises as an empirical measure of the form , where Y _j are independent with distribution Λ _i. More generally, can be a normalised point process with mean measure τΛ _i, i.e.

This encapsulates the case of empirical measure when τ is an integer and is a binomial point process. The parameter τ is the expected number of observed points over the entire space ; the larger τ is, the more information gives on Λ _i.

Except if is an empirical measure, there is one difficulty in the above setting that needs to be addressed. Unless is binomial, there is a positive probability that and no points pertaining to Λ _i are observed. In the asymptotic setup below, conditions will be imposed to ensure that this probability becomes negligible as n →∞. For concreteness we define for some fixed measure λ ⁽⁰⁾ that will be of minor importance. This can be a Dirac measure at 0, a certain fixed Gaussian measure, or (normalised) Lebesgue measure on some bounded set in case . We can now replace the estimator λ _n by , defined as any minimiser of

which exists by Corollary 3.​1.​3.

As a generalisation of the discrete case discussed in Sect. 1.​3, the Fréchet mean of discrete measures can be computed exactly. Suppose that is nonzero for all i. Then each is a discrete measure supported on N _i points. One can then recast the multimarginal formulation (see Sect. 3.​1.​2) as a finite linear program, solve it, and “average” the solution as in Proposition 3.​1.​2 in order to obtain (an alternative linear programming formulation for finding a Fréchet mean is given by Anderes et al. [14]). Thus, can be computed in finite time, even when is infinite-dimensional.

Finally, a remark about measurability is in order. Point processes can be viewed as random elements in endowed with the vague topology induced from convergence of integrals of continuous functions with compact support. If μ _n converge to μ vaguely, and a _n are numbers that converge to a, then a _nμ _n → aμ vaguely. Thus, is a continuous function of the pair and can be viewed as a random measure with respect to the vague topology. The restriction of the vague topology to probability measures is equivalent to the weak topology,⁵ and therefore vague, weak, and Wasserstein measurability are all equivalent.

4.3.3 Smoothing

Even when the computational complexity involved in calculating is tractable, there is another reason not to use it as an estimator for λ. If one has a-priori knowledge that λ is smooth, it is often desirable to estimate it by a smooth measure. One way to achieve this would be to apply some smoothing technique to using, e.g., kernel density estimation. However, unless the number of observed points from each measure is the same N ₁ = ⋯ = N _n = N, will usually be concentrated on many points, essentially N ₁ + ⋯ + N _n of them. In other words, the Fréchet mean is concentrated on many more points than each of the measures , thus potentially hindering its usefulness as a mean because it will not be a representative of the sample.

This is most easily seen when , in which case each is a discrete uniform measure on points , where we assume for simplicity that the points are not repeated (this will happen with probability one if Λ _i is diffuse). If we now set G _i to be the distribution function of , then the quantile function is piecewise constant on each interval (k, k + 1]∕N _i with jumps at

The Fréchet mean has quantile function and will have jumps at every point of the form k∕N _i for k ≤ N _i and i = 1, …, n. In the worst-case scenario, when no pair from N _i has a common divisor, there will be

jumps for G ⁻¹, which is the number of points on which the Fréchet mean will be supported. (All the ’s have a jump at one which thus needs to be counted once rather than n times.)

By counting the number of redundancies in the constraints matrix of the linear program, one can show that this is in general an upper bound on the number of support points of the Fréchet mean.

An alternative approach is to first smooth each observation and then calculate the Fréchet mean. Since it is easy to bound the Wasserstein distances when dealing with convolutions, we will employ kernel density estimation, although other smoothing approaches could be used as well.

To simplify the exposition, we provide the technical details only when , but a similar construction will work when the dimension of is infinite. Let be a continuous, bounded, strictly positive isotropic density function with unit variance: ψ(x) = ψ ₁(∥x∥) with ψ ₁ nonincreasing and

(Besides the boundedness all these properties can be relaxed, and if even boundedness is not necessary.) A classical example for ψ is the standard Gaussian density in . Define the rescaled version ψ _σ(x) = σ ^−dψ(x∕σ) for all σ > 0. We can then replace by a smooth proxy . If is a sum of Dirac masses at , then

If N _i = 0, one can either use λ ⁽⁰⁾ or λ ⁽⁰⁾ ∗ ψ _σ; this event will have negligible probability anyway.

For the purpose of approximating , this convolution is an acceptable estimator, because as was seen in the proof of Theorem 2.​2.​7,

But the measure has a strictly positive density throughout . If we know that Λ is supported on a convex compact , it is desirable to construct an estimator that has the same support K. The first idea that comes to mind is to project to K (see Proposition 3.​2.​4), as this will further decrease the Wasserstein distance, but the resulting measure will then have positive mass on the boundary of K, and will not be absolutely continuous. We will therefore use a different strategy: eliminate all the mass outside K and redistribute it on K. The simplest way to do this is to restrict to K and renormalise the restriction to be a probability measure. For technical reasons, it will be more convenient to bound the Wasserstein distance when the restriction and renormalisation is done separately on each point of . This yields the measure

(4.2)

Lemma 4.4.2 below shows that for some finite constant C. It is apparent that is a continuous function of and σ, so is measurable; in any case this is not particularly important because σ will vanish, so asymptotically and the latter is measurable.

Our final estimator for λ is defined as the minimiser of

Since the measures are absolutely continuous, is unique. We refer to as the regularised Fréchet–Wasserstein estimator, where the regularisation comes from the smoothing and the possible restriction to K.

In the case , can be constructed via averaging of quantile functions. Let be the distribution function of . Then is the measure with quantile function

and distribution function

By construction, the are continuous and strictly increasing, so the inverses are proper inverses and one does not to use the right-continuous inverse as in Sect. 3.​1.​4.

If and d ≥ 2, then there is no explicit expression for , although it exists and is unique. In the next chapter, we present a steepest descent algorithm that approximately constructs by taking advantage of the differentiability properties of the Fréchet functional in Sect. 3.​1.​6.

4.3.4 Estimation of Warpings and Registration Maps

Once estimators , i = 1, …, n and are constructed, it is natural to estimate the map and its inverse (when Λ _i are absolutely continuous; see the discussion after Assumptions 3 below) by the plug-in estimators

The latter, the registration maps, can then be used in order to register the points Π _i via

It is thus reasonable to expect that if is a good estimator, then its composition with T _i should be close to the identity and should be close to Π _i.

4.3.5 Unbiased Estimation When

In the same way, Fréchet means extend the notion of mean to non-Hilbertian spaces, they also extend the definition of unbiased estimators. Let H be a separable Hilbert space (or a convex subset thereof) and suppose that is a random element in H whose distribution μ _θ depends on a parameter θ ∈ H. Then is unbiased for θ if for all θ ∈ H

(We use the standard notation in the sequel.) This is equivalent to

In view of that, one can define unbiased estimators of as measurable functions δ = δ(Λ ₁, …, Λ _n) for which

This definition was introduced by Lehmann [88].

Unbiased estimators allow us to avoid the problem of over-registering (the so-called pinching effect; Kneip and Ramsay [82, Section 2.4]; Marron et al. [90, p. 476]). An extreme example of over-registration is if one “aligns” all the observed patterns into a single fixed point x ₀. The registration will then seem “successful” in the sense of having no residual phase variation, but the estimation is clearly biased because the points are not registered to the correct reference measure. Thus, requiring the estimator to be unbiased is an alternative to penalising the registration maps.

Due to the Hilbert space embedding of , it is possible to characterise unbiased estimators in terms of a simple condition on their quantile functions. As a corollary, λ _n, the Fréchet mean of {Λ ₁, …, Λ _n}, is unbiased. Our regularised Fréchet–Wasserstein estimator can then be interpreted as approximately unbiased, since it approximates the unobservable λ _n.

Proposition 4.3.1 (Unbiased Estimators in )

Let Λ be a random measure in with finite Fréchet functional and let λ be the unique Fréchet mean of Λ (Theorem 3.2.11 ). An estimator δ constructed as a function of a sample (Λ ₁, …, Λ _n) is unbiased for λ if and only if the left-continuous representatives (in L ₂(0, 1)) satisfy for all x ∈ (0, 1).

Proof

The proof is straightforward from the definition: δ is unbiased if and only if for all λ and all γ,

which is equivalent to . In other words, these two functions must equal almost everywhere on (0, 1), and their left-continuous representatives must equal everywhere (the fact that has such a representative was established in Sect. 3.​1.​4).

To show that δ = λ _n is unbiased, we simply invoke Theorem 3.​2.​11 twice to see that

which proves unbiasedness of δ.

4.4 Consistency

In functional data analysis, one often assumes that the number of curves n and the number of observed points per curve m both diverge to infinity. An analogous framework for point processes would similarly require the number of point processes n as well as the expected number of points τ per processes to diverge. A technical complication arises, however, because the mean measures do not suffice to characterise the distribution of the processes. Indeed, if one is given a point processes Π with mean measure λ (not necessarily a probability measure), and τ is an integer, there is no unique way to define a process Π ^(τ) with mean measure τλ. One can define Π ^(τ) = τΠ, so that every point in Π will be counted τ times. Such a construction, however, can never yield a consistent estimator of λ, even when τ →∞.

Another way to generate a point process with mean measure τλ is to take a superposition of τ independent copies of Π. In symbols, this means

with (Π _i) independent, each having the same distribution as Π. This superposition scheme gives the possibility to use the law of large numbers. If τ is not an integer, then this construction is not well-defined but can be made so by assuming that the distribution of Π is infinitely divisible. The reader willing to assume that τ is always an integer can safely skip to Sect. 4.4.1; all the main ideas are developed first for integer values of τ and then extended to the general case.

A point process Π is infinitely divisible if for every integer m there exists a collection of m independent and identically distributed such that

If Π is infinitely divisible and τ = k∕m is rational, then can define π ^(τ) using km independent copies of Π ^(1∕m):

One then deals with irrational τ via duality and continuity arguments, as follows. Define the Laplace functional of Π by

The Laplace functional characterises the distribution of the point process, generalising the notion of Laplace transform of a random variable or vector (Karr [79, Theorem 1.12]). By definition, it translates convolutions into products. When Π = Π ⁽¹⁾ is infinitely divisible, the Laplace functional L ₁ of Π takes the form (Kallenberg [75, Chapter 6]; Karr [79, Theorem 1.43])

The Laplace functional of Π ^(τ) is L _τ(f) = [L ₁(f)]^τ for any rational τ, which simply amounts to multiplying the measure ρ by the scalar τ. One can then do the same for an irrational τ, and the resulting Laplace functional determines the distribution of Π ^(τ) for all τ > 0.

4.4.1 Consistent Estimation of Fréchet Means

We are now ready to define our asymptotic setup. The following assumptions will be made. Notice that the Wasserstein geometry does not appear explicitly in these assumptions, but is rather derived from them in view of Theorem 4.2.4. The compactness requirement can be relaxed under further moment conditions on λ and the point process Π; we focus on the compact case for the simplicity and because in practice the point patterns will be observed on a bounded observation window.

Assumptions 3

Let be a compact convex nonempty set, λ an absolutely continuous probability measure on K, and τ _n a sequence of positive numbers. Let Π be a point processes on K with mean measure λ. Finally, define U = intK.

• For every n, let be independent point processes, each having the same distribution as a superposition of τ _n copies of Π.

• Let T be a random injective function on K (viewed as a random element in C _b(K, K) endowed with the supremum norm) such that T(x) ∈ U for x ∈ U (that is, T ∈ C _b(U, U)) with nonsingular derivative for almost all x ∈ U, that is a gradient of a convex function. Let {T ₁, …, T _n} be independent and identically distributed as T.

• For every x ∈ U, assume that .

• Assume that the collections and are independent.

• Let be the warped point processes, having conditional mean measures .

• Define by the smoothing procedure (4.2), using bandwidth (possibly random).

The dependence of the estimators on n will sometimes be tacit. But Λ _idoes not depend on n.

By virtue of Theorem 4.2.4, λ is a Fréchet mean of the random measure Λ = T#λ. Uniqueness of this Fréchet mean will follow from Proposition 3.​2.​7 if we show that Λ is absolutely continuous with positive probability. This is indeed the case, since T is injective and has a nonsingular Jacobian matrix; see Ambrosio et al. [12, Lemma 5.5.3]. The Jacobian assumption can be removed when , because Fréchet means are always unique by Theorem 3.​2.​11.

Notice that there is no assumption about the dependence between rows. Assumptions 3 thus cover, in particular, two different scenarios:

• Full independence: here the point processes are independent across rows, that is, and are also independent.

• Nested observations: here includes the same points as and additional points, that is, is a superposition of and another point process distributed as (τ _n+1 − τ _n)Π.

Needless to say, Assumptions 3 also encompass binomial processes when τ _n are integers, as well as Poisson processes or, more generally, Poisson cluster processes.

We now state and prove the consistency result for the estimators of the conditional mean measures Λ _i and the structural mean measure λ.

Theorem 4.4.1 (Consistency)

If Assumptions 3 hold, almost surely and τ _n →∞ as n →∞, then:

1. 1.
   The estimators defined by (4.2), constructed with bandwidth , are Wasserstein-consistent for the conditional mean measures: for all i such that

   
    
2. 2.
   The regularised Fréchet–Wasserstein estimator of the structural mean measure (as described in Sect. 4.3) is strongly Wasserstein-consistent,

   
    

Convergence in 1. holds almost surely under the additional conditions that and . If σ _n → 0 only in probability, then convergence in 2. still holds in probability.

Theorem 4.4.1 still holds without smoothing (σ _n = 0). In that case, is possibly not unique, and the theorem should be interpreted in a set-valued sense (as in Proposition 1.​7.​8): almost surely, any choice of minimisers converges to λ as n →∞.

The preceding paragraph notwithstanding, we will usually assume that some smoothing is present, in which case is unique and absolutely continuous by Proposition 3.​1.​8. The uniform Lipschitz bounds for the objective function show that if we restrict the relevant measures to be absolutely continuous, then is a continuous function of and hence is measurable; this is again a minor issue because many arguments in the proof hold for each ω ∈ Ω separately. Thus, even if is not measurable, the proof shows that the convergence holds outer almost surely or in outer probability.

The first step in proving consistency is to show that the Wasserstein distance between the unsmoothed and the smoothed estimators of Λ _i vanishes with the smoothing parameter. The exact rate of decay will be important to later establish the rate of convergence of to λ, and is determined next.

Lemma 4.4.2 (Smoothing Error)

There exists a finite constant C _ψ,K, depending only on ψ and on K, such that

(4.3)

Since the smoothing parameter will anyway vanish, this restriction to small values of σ is not binding. The constant C _ψ,K is explicit. When , a more refined construction allows to improve this constant in some situations, see Panaretos and Zemel [100, Lemma 1].

Proof

The idea is that (4.2) is a sum of measures with mass 1∕N _i that can be all sent to the relevant point x _j, and we refer to page 98 in the supplement for the precise details.

Proof(Proof of Theorem 4.4.1) The proof, detailed on page 97 of the supplement, follows the following steps: firstly, one shows the convergence in probability of to Λ _i. This is basically a corollary of Karr [79, Proposition 4.8] and the smoothing bound from Lemma 4.4.2.

To prove claim (2) one considers the functionals, defined on :

Since K is compact, they are all locally Lipschitz, so their differences can be controlled by the distances between Λ _i, , and . The first distance vanishes since the intensity τ →∞, and the second by the smoothing bound. Another compactness argument yields that uniformly on , and so the minimisers converge.

The almost sure convergence in (1) is proven as follows. Under the stronger conditions at the end of the theorem’s statement, for any fixed ,

by the law of large numbers. This extends to all rational a’s, then to all a by approximation. The smoothing error is again controlled by Lemma 4.4.2.

4.4.2 Consistency of Warp Functions and Inverses

We next discuss the consistency of the warp and registration function estimators. These are key elements in order to align the observed point patterns . Recall that we have consistent estimators for Λ _i and for λ. Then is estimated by and is estimated by . We will make the following extra assumptions that lead to more transparent statements (otherwise one needs to replace K with the set of Lebesgue points of the supports of λ and Λ _i).

Assumptions 4 (Strictly Positive Measures)

In addition to Assumptions 3suppose that:

1. 1.

   λ has a positive density on K (in particular, suppλ = K);

    
2. 2.

   T is almost surely surjective on U = intK (thus a homeomorphism of U).

    

As a consequence almost surely.

Theorem 4.4.3 (Consistency of Optimal Maps)

Let Assumptions 4 be satisfied in addition to the hypotheses of Theorem 4.4.1. Then for any i such that and any compact set S ⊆intK,

Almost sure convergence can be obtained under the same provisions made at the end of the statement of Theorem 4.4.1.

A few technical remarks are in order. First and foremost, it is not clear that the two suprema are measurable. Even though T _i and are random elements in , their estimators are only defined in an L ₂ sense. The proof of Theorem 4.4.3 is done ω-wise. That is, for any ω in the probability space such that Theorem 4.4.1 holds, the two suprema vanish as n →∞. In other words, the convergence holds in outer probability or outer almost surely.

Secondly, assuming positive smoothing, the random measures are smooth with densities bounded below on K, so are defined on the whole of U (possibly as set-valued functions on a Λ _i-null set). But the only known regularity result for is an upper bound on its density (Proposition 3.​1.​8), so it is unclear what is its support and consequently what is the domain of definition of .

Lastly, when the smoothing parameter σ is zero, and are not defined. Nevertheless, Theorem 4.4.3 still holds in the set-valued formulation of Proposition 1.​7.​11, of which it is a rather simple corollary:

Proof(Proof of Theorem 4.4.3) The proof amounts to setting the scene in order to apply Proposition 1.​7.​11 of stability of optimal maps. We define

and verify the conditions of the proposition. The weak convergence of μ _n to μ and ν _n to ν is the conclusion of Theorem 4.4.1; the finiteness is apparent because K is compact and the uniqueness follows from the assumed absolute continuity of Λ _i. Since in addition is uniquely defined on U = intK which is an open convex set, the restrictions on Ω in Proposition 1.​7.​11 are redundant. Uniform convergence of to T _i is proven in the same way.

Corollary 4.4.4 (Consistency of Point Pattern Registration)

For any i such that ,

The division by the number of observed points ensures that the resulting measures are probability measures; the relevant information is contained in the point patterns themselves, and is invariant under this normalisation.

Proof

The law of large numbers entails that , so in particular is almost surely not zero when n is large. Since , we have the upper bound

Fix a compact Ω ⊆intK and split the integral to Ω and its complement. Then

by the law of large numbers, where d _K is the diameter of K. By writing intK as a countable union of compact sets (and since λ is absolutely continuous), this can be made arbitrarily small by choice of Ω.

We can easily bound the integral on Ω by

But T _i(Ω) is a compact subset of U = intK, because T _i ∈ C _b(U, U). The right-hand side therefore vanishes as n →∞ by Theorem 4.4.3, and this completes the proof.

Possible extensions pertaining to the boundary of K are discussed on page 33 of the supplement.

4.5 Illustrative Examples

In this section, we illustrate the estimation framework put forth in this chapter by considering an example of a structural mean λ with a bimodal density on the real line. The unwarped point patterns Π originate from Poisson processes with mean measure λ and, consequently, the warped points are Cox processes (see Sect. 4.1.2). Another scenario involving triangular densities can be found in Panaretos and Zemel [100].

4.5.1 Explicit Classes of Warp Maps

As a first step, we introduce a class of random warp maps satisfying Assumptions 2, that is, increasing maps that have as mean the identity function. The construction is a mixture version of similar maps considered by Wang and Gasser [128, 129].

For any integer k define ζ _k : [0, 1] → [0, 1] by

(4.4)

Clearly ζ _k(0) = 0, ζ _k(1) = 1 and ζ _k is smooth and strictly increasing for all k. Figure 4.4a plots ζ _k for k = −3, …, 3. To make ζ _k a random function, we let k be an integer-valued random variable. If the latter is symmetric, then we have

By means of mixtures, we replace this discrete family by a continuous one: let J > 1 be an integer and V = (V ₁, …, V _J) be a random vector following the flat Dirichlet distribution (uniform on the set of nonnegative vectors with v ₁ + ⋯ + v _J = 1). Take independently k _j following the same distribution as k and define

(4.5)

Since V _j is positive, T is increasing and as (V _j) sums up to unity T has mean identity. Realisations of these warp functions are given in Fig. 4.4b and c for J = 2 and J = 10, respectively. The parameters (k _j) were chosen as symmetrised Poisson random variables: each k _j has the law of XY with X Poisson with mean 3 and for Y and X independent. When J = 10 is large, the function T deviates only mildly from the identity, since a law of large numbers begins to take effect. In contrast, J = 2 yields functions that are quite different from the identity. Thus, it can be said that the parameter J controls the variance of the random warp function T.

Fig. 4.4

(a) The functions {ζ ₋₃, …, ζ ₃}; (b) realisations of T defined by (4.5) with J = 2 and k _j symmetrisations of Poisson random variables with mean 3; (c) realisations of T defined by (4.5) with J = 10 and k _j as in (b)

4.5.2 Bimodal Cox Processes

Let the structural mean measure λ be a mixture of a bimodal Gaussian distribution (restricted to K = [−16, 16]) and a beta background on the interval [−12, 12], so that mass is added at the centre of K but not near the boundary. In symbols this is given as follows. Let φ be the standard Gaussian density and let β _α,β denote the density of a the beta distribution with parameters α and β. Then λ is chosen as the measure with density

(4.6)

where 𝜖 ∈ [0, 1] is the weight of the beta background. (We ignore the loss of a negligible amount of mass due to the restriction of the Gaussians to [−16, 16].) Plots of the density and distribution functions are given in Fig. 4.5.

Fig. 4.5

Density and distribution functions corresponding to (4.6) with 𝜖 = 0 and 𝜖 = 0.15

The main criterion for the quality of our regularised Fréchet–Wasserstein estimator will be its success in discerning the two modes at ± 8; these will be smeared by the phase variation arising from the warp functions.

We next simulated 30 independent Poisson processes with mean measure λ, 𝜖 = 0.1, and total intensity (expected number of points) τ = 93. In addition, we generated warp functions as in (4.5) but rescaled to [−16, 16]; that is, having the same law as the functions

from K to K. These cause rather violent phase variation, as can be seen by the plots of the densities and distribution functions of the conditional measures Λ = T#λ presented in Fig. 4.6a and b; the warped points themselves are displayed in Fig. 4.6c.

Fig. 4.6

(a) 30 warped bimodal densities, with density of λ given by (4.6) in solid black; (b) their corresponding distribution functions, with that of λ in solid black; (c) 30 Cox processes, constructed as warped versions of Poisson processes with mean intensity 93f using as warp functions the rescaling to [−16.16] of (4.5)

Using these warped point patterns, we construct the regularised Fréchet– Wasserstein estimator employing the procedure described in Sect. 4.3. Each was smoothed with a Gaussian kernel and bandwidth chosen by unbiased cross validation. We deviate slightly from the recipe presented in Sect. 4.3 by not restricting the resulting estimates to the interval [−16, 16], but this has no essential effect on the finite sample performance. The regularised Fréchet–Wasserstein estimator serves as the estimator of the structural mean λ and is shown in Fig. 4.7a. It is contrasted with λ at the level of distribution functions, as well as with the empirical arithmetic mean; the latter, the naive estimator, is calculated by ignoring the warping and simply averaging linearly the (smoothed) empirical distribution functions across the observations. We notice that is rather successful at locating the two modes of λ, in contrast with the naive estimator that is more diffuse. In fact, its distribution function increases approximately linearly, suggesting a nearly constant density instead of the correct bimodal one.

Fig. 4.7

(a) Comparison between the regularised Fréchet–Wasserstein estimator, the empirical arithmetic mean, and the true distribution function, including residual curves centred at y = 3∕4; (b) The estimated warp functions; (c) Kernel estimates of the density function f of the structural mean, based on the warped and registered point patterns

Estimators of the warp maps , depicted in Fig. 4.7b, and their inverses, are defined as the optimal maps between and the estimated conditional mean measures, as explained in Sect. 4.3.4. Then we register the point patterns by applying to them the inverse estimators (Fig. 4.8). Figure 4.7c gives two kernel estimators of the density of λ constructed from a superposition of all the warped points and all the registered ones. Notice that the estimator that uses the registered points is much more successful than the one using the warped ones in discerning the two density peaks. This is not surprising after a brief look at Fig. 4.8, where the unwarped, warped, and registered points are displayed. Indeed, there is very high concentration of registered points around the true location of the peaks, ± 8. This is not the case for the warped points because of the phase variation that translates the centres of concentration for each individual observation. It is important to remark that the fluctuations in the density estimator in Fig. 4.7c are not related to the registration procedure, and could be reduced by a better choice of bandwidth. Indeed, our procedure does not attempt to estimate the density, but rather the distribution function.

Fig. 4.8

Bimodal Cox processes: (a) the observed warped point processes; (b) the unobserved original point processes; (c) the registered point processes

Figure 4.9 presents a superposition of the regularised Fréchet–Wasserstein estimators for 20 independent replications of the experiment, contrasted with a similar superposition for the naive estimator. The latter is clearly seen to be biased around the two peaks, while the regularised Fréchet–Wasserstein seems approximately unbiased, despite presenting fluctuations. It always captures the bimodal nature of the density, as is seen from the two clear elbows in each realisation.

Fig. 4.9

(a) Sampling variation of the regularised Fréchet–Wasserstein mean and the true mean measure λ for 20 independent replications of the experiment; (b) sampling variation of the arithmetic mean, and the true mean measure λ for the same 20 replications; (c) superposition of (a) and (b). For ease of comparison, all three panels include residual curves centred at y = 3∕4

To illustrate the consistency of the regularised Fréchet–Wasserstein estimator for λ as shown in Theorem 4.4.1, we let the number of processes n as well as the expected number of observed point per process τ to vary. Figures 4.10 and 4.11 show the sampling variation of for different values of n and τ. We observe that as either of these increases, the realisations indeed approach λ. The figures suggest that, in this scenario, the amplitude variation plays a stronger role than the phase variation, as the effect of τ is more substantial.

Fig. 4.10

Sampling variation of the regularised Fréchet–Wasserstein mean and the true mean measure λ for 20 independent replications of the experiment, with 𝜖 = 0 and n = 30. Left: τ = 43; middle: τ = 93; right: τ = 143. For ease of comparison, all three panels include residual curves centred at y = 3∕4

Fig. 4.11

Sampling variation of the regularised Fréchet–Wasserstein mean and the true mean measure λ for 20 independent replications of the experiment, with 𝜖 = 0 and τ = 93. Left: n = 30; middle: n = 50; right: n = 70. For ease of comparison, all three panels include residual curves centred at y = 3∕4.

4.5.3 Effect of the Smoothing Parameter

In order to work with measures of strictly positive density, the observed point patterns have been smoothed using a kernel function. This necessarily incurs an additional bias that depends on the bandwidth σ _i. The asymptotics (Theorem 4.4.1) guarantee the consistency of the estimators, in particular the regularised Fréchet–Wasserstein estimator , provided that . In our simulations, we choose σ _i in a data-driven way by employing unbiased cross validation. To gauge for the effect of the smoothing, we carry out the same estimation procedure but with σ _i multiplied by a parameter s. Figure 4.12 presents the distribution function of as a function of s. Interestingly, the curves are nearly identical as long as s ≤ 1, whereas when s > 1, the bias becomes more substantial.

Fig. 4.12

Regularised Fréchet–Wasserstein mean as a function of the smoothing parameter multiplier s, including residual curves. Here, n = 30 and τ = 143

These findings are reaffirmed in Fig. 4.13 that show the registered point processes again as a function of s. We see that only minor differences are present as s varies from 0.1 to 1, for example, in the grey (8), black (17), and green (19) processes. When s = 3, the distortion becomes quite more substantial. This phenomenon repeats itself across all combinations of n, τ, and s tested.

Fig. 4.13

Registered point processes as a function of the smoothing parameter multiplier s. Left: s = 0.1; middle: s = 1; right: s = 3. Here, n = 30 and τ = 43

4.6 Convergence Rates and a Central Limit Theorem on the Real Line

Since the conditional mean measures Λ _i are discretely observed, the rate of convergence of our estimators will be affected by the rate at which the number of observed points per process increases to infinity. The latter is controlled by the next lemma, which is valid for any complete separable metric space .

Lemma 4.6.1 (Number of Points Grows Linearly)

Let denote the total number of observed points. If , then there exists a constant C _Π > 0, depending only on the distribution of Π, such that almost surely

In particular, there are no empty point processes, so the normalisation is well-defined. If Π is a Poisson process, then we have the more precise result

Remark 4.6.2

One can also show that the limit superior of the same quantity is bounded by a constant . If is bounded below, then the same result holds but with worse constants. If only τ _n →∞, then the result holds for each i separately but in probability.

The proof is a simple application of Chernoff bounds; see page 108 in the supplement.

With Lemma 4.6.1 under our belt, we can replace terms of the order by the more transparent order τ _n. As in the consistency proof, the idea is to write

and control each term separately. The first term corresponds to the phase variation, and comes from the approximation of the theoretical expectation F by a sample mean F _n. The second term is associated with the amplitude variation resulting from observing Λ _i discretely. The third term can be viewed as the bias incurred by the smoothing procedure. Accordingly, the rate at which converges to λ is a sum of three separate terms. We recall the standard terminology: if X _n and Y _n are random variables, then means that the sequence (X _n∕Y _n) is bounded in probability, which by definition is the condition

Instead of , we will sometimes write . The former notation emphasises the condition that X _n grows no faster than Y _n, while the latter stresses that Y _n grows at least as fast as X _n (which is of course the same assertion). Finally, means that X _n∕Y _n → 0 in probability.

Theorem 4.6.3 (Convergence Rates on )

Suppose in addition to Assumptions 3 that d = 1, and that Π is either a Poisson process or a binomial process. Then

where all the constants in the terms are explicit.

Remark 4.6.4

Unlike classical density estimation, no assumptions on the rate of decay of σ _n are required, because we only need to estimate the distribution function and not the derivative. If the smoothing parameter is chosen to be for some α > 0 and , then by Lemma 4.6.1 . For example, if Rosenblatt’s rule α = 1∕5 is employed, then the term can be replaced by .

One can think about the parameter τ as separating the sparse and dense regimes as in classical functional data analysis (see also Wu et al. [132]). If τ is bounded, then the setting is ultra sparse and consistency cannot be achieved. A sparse regime can be defined as the case where τ _n →∞ but slower than . In that case, consistency is guaranteed, but some point patterns will be empty. The dense regime can be defined as τ _n ≫ n ², in which case the amplitude variation is negligible asymptotically when compared with the phase variation.

The exponent − 1∕4 of τ _n can be shown to be optimal without further assumptions, but it can be improved to − 1∕2 if for some 𝜖 > 0, where f _Λ is the density of Λ (see Sect. 4.7). In terms of T, the condition is that for some 𝜖 and λ has a density bounded below. When this is the case, τ _n needs to compared with n rather than n ² in the next paragraph and the next theorem.

Theorem 4.6.3 provides conditions for the optimal parametric rate to be achieved: this happens if we set σ _n to be of the order or less and if τ _n is of the order n ² or more. But if the last two terms in Theorem 4.6.3 are negligible with respect to n ^−1∕2, then a sort of central limit theorem holds for :

Theorem 4.6.5 (Asymptotic Normality)

In addition to the conditions of Theorem 4.6.3, assume that τ _n∕n ² →∞, and λ possesses an invertible distribution function F _λ on K. Then

for a zero-mean Gaussian process Z with the same covariance operator of T (the latter viewed as a random element in L ₂(λ)), namely with covariance kernel

If the density f _λ exists and is (piecewise) continuous and bounded below on K, then the weak convergence also holds in L ₂(K).

In view of Sect. 2.​3, Theorem 4.6.5 can be interpreted as asymptotic normality of in the tangential sense: converges to a Gaussian random element in the tangent space Tan_λ, which is a subset of L ₂(λ). The additional smoothness conditions allow to switch to the space L ₂(K), which is independent of the unknown template measure λ.

See pages 109 and 110 in the supplement for detailed proofs of these theorems. Below we sketch the main ideas only.

Proof(Proof of Theorem 4.6.3)

The quantile formula from Sect. 1.​5 and the average quantile formula for the Fréchet mean (Sect. 3.​1.​4) show that the oracle empirical mean follows a central limit theorem in L ₂(0, 1). Since we work in the Hilbert space L ₂(0, 1), Fréchet means are simple averages, so the errors in the Fréchet mean have the same rate as the errors in the Fréchet functionals. The smoothing term is easily controlled by Lemma 4.4.2.

Controlling the amplitude term is more difficult. Bounds can be given using the machinery sketched in Sect. 4.7, but we give a more elementary proof by reducing to the 1-Wasserstein case (using (2.​2)), which can be more easily handled in terms of distributions functions (Corollary 1.​5.​3).

Proof(Proof of Theorem 4.6.5) The hypotheses guarantee that the amplitude and smoothing errors are negligible and

where GP is the Gaussian process defined in the proof of Theorem 4.6.3. One then employs a composition with F _λ.

4.7 Convergence of the Empirical Measure and Optimality

One may find the term in Theorem 4.6.3 to be somewhat surprising, and expect that it ought to be . The goal of this section is to show why the rate is optimal without further assumptions and discuss conditions under which it can be improved to the optimal rate . For simplicity, we concentrate on the case τ _n = n and assume that the point process Π is binomial; the Poisson case being easily obtained from the simplified one (using Lemma 4.6.1). We are thus led to study rates of convergence of empirical measures in the Wasserstein space. That is to say, for a fixed exponent p ≥ 1 and a fixed measure , we consider independent random variables X ₁, … with law μ and the empirical measure . The first observation is that :

Lemma 4.7.1

Let be any measure. Then

Proof

This result has been established in an almost sure sense in Proposition 2.​2.​6. To extend to convergence in expectation observe that

Thus, the random variable is bounded by the sample average Z _n of a random variable that has a finite expectation. A version of the dominated converge theorem (given on page 111 in the supplement) implies that . Now invoke Jensen’s inequality.

Remark 4.7.2

The sequence is not monotone, as the simple example μ = (δ ₀ + δ ₁)∕2 shows (see page 111 in the supplement).

The next question is how quickly vanishes when . We shall begin with two simple general lower bounds, then discuss upper bounds in the one-dimensional case, put them in the context of Theorem 4.6.3, and finally briefly touch the d-dimensional case.

Lemma 4.7.3 ( Lower Bound)

Let be nondegenerate. Then there exists a constant c(μ) > 0 such that for all p ≥ 1 and all n

Proof

Let X ∼ μ and let a≠b be two points in the support μ. Consider , a bounded 1-Lipschitz function such that f(a) = 0 < f(b). Then

by the central limit theorem and the Kantorovich–Rubinstein theorem (1.​11).

For discrete measures, the rates scale badly with p. More generally:

Lemma 4.7.4 (Separated Support)

Suppose that there exist Borel sets such that μ(A ∪ B) = 1,

Then for any p ≥ 1 there exists c _p(μ) > 0 such that .

Any nondegenerate finitely discrete measure μ satisfies this condition, and so do “non-pathological” countably discrete ones. (An example of a “pathological” measure is one assigning positive mass to any rational number.)

Proof

Let k ∼ B(n, q = μ(A)) denote the number of points from the sample (X ₁, …, X _n) that fall in A. Then a mass of |k∕n − q| must travel between A and B, a distance of at least . Thus, , and the result follows from the central limit theorem for k; see page 112 in the supplement for the full details.

These lower bounds are valid on any separable metric space. On the real line, it is easy to obtain a sufficient condition for the optimal rate n ^−1∕2 to be attained for W ₁: since F _n(t) ∼ B(n, F(t)) has variance F(t)(1 − F(t))∕n, we have (by Fubini’s theorem and Jensen’s inequality)

so that W ₁(μ _n, μ) is of the optimal order n ^−1∕2 if

Since the integrand is bounded by 1∕2, this is certainly satisfied if μ is compactly supported. The J ₁ condition is essentially a moment condition, since for any δ > 0, we have for X ∼ μ that . It turns out that this condition is necessary, and has a more subtle counterpart for any p ≥ 1. Let f denote the density of the absolutely continuous part of μ (so f ≡ 0 if μ is discrete).

Theorem 4.7.5 (Rate of Convergence of Empirical Measures)

Let p ≥ 1 and . The condition

is necessary and sufficient for .

See Bobkov and Ledoux [25, Theorem 5.10] for a proof for the J _p condition, and Theorems 5.1 and 5.3 for the values of the constants and a stronger result.

When p > 1, for J _p(μ) to be finite, the support of μ must be connected; this is not needed when p = 1. Moreover, the J _p condition is satisfied when f is bounded below (in which case the support of μ must be compact). However, smoothness alone does not suffice, even for measures with positive density on a compact support. More precisely, we have:

Proposition 4.7.6

For any rate 𝜖 _n → 0 there exists a measure μ on [−1, 1] with positive C ^∞ density there, and such that for all n

The rate n ^−1∕(2p) from Lemma 4.7.4 is the worst among compactly supported measures on . Indeed, by Jensen’s inequality and (2.​2), for any μ ∈ P([0, 1]),

The proof of Proposition 4.7.6 is done by “smoothing” the construction in Lemma 4.7.4, and is given on page 113 in the supplement.

Let us now put this in the context of Theorem 4.6.3. In the binomial case, since each and each Λ _i are independent, we have

(In the Poisson case, we need to condition on and then estimate its inverse square root as is done in the proof of Theorem 4.6.3.) Therefore, a sufficient condition for the rate to hold is that and a necessary condition is that . These hold if there exists δ > 0 such that with probability one Λ has a density bounded below by δ. Since Λ = T#λ, this will happen provided that λ itself has a bounded below density and T has a bounded below derivative. Bigot et al. [23] show that the rate cannot be improved.

We conclude by proving a lower bound for absolutely continuous measures and stating, without proof, an upper bound.

Proposition 4.7.7

Let have an absolutely continuous part with respect to Lebesgue measure, and let ν _n be any discrete measure supported on n points (or less). Then there exists a constant C(μ) > 0 such that

Proof

Let f be the density of the absolutely continuous part μ _c, and observe that for some finite number M,

Let x ₁, …, x _n be the support points of ν _n and 𝜖 > 0. Let μ _c,M be the restriction of μ _c to the set where the density is smaller than M. The union of balls B _𝜖(x _i) has μ _c,M-measure of at most

if 𝜖 ^d = δ(nMC _d)⁻¹. Thus, a mass 2δ − δ = δ must travel more than 𝜖 from ν _n to μ in order to cover μ _c,M. Hence

The lower bound holds because we need 𝜖 ^−d balls of radius 𝜖 in order to cover a sufficiently large fraction of the mass of μ. The determining quantity for upper bounds on the empirical Wasserstein distance is the covering numbers

Since μ is tight, these are finite for all 𝜖, τ > 0, and they increase as 𝜖 and τ approach zero. To put the following bound in context, notice that if μ is compactly supported on , then N(μ, 𝜖, 0) ≤ K𝜖 ^−d.

Theorem 4.7.8

If for some d>2p, N(μ, 𝜖, 𝜖 ^{dp∕(d−2p)})≤K𝜖 ^−d, then .

Comparing this with the lower bound in Lemma 4.7.4, we see that in the high-dimensional regime d > 2p, absolutely continuous measures have a worse rate than discrete ones. In the low-dimensional regime d < 2p, the situation is opposite. We also obtain that for d > 2 and a compactly supported absolutely continuous , .

4.8 Bibliographical Notes

Our exposition in this chapter closely follows the papers Panaretos and Zemel [100] and Zemel and Panaretos [134].

Books on functional data analysis include Ramsay and Silverman [109, 110], Ferraty and Vieu [51], Horváth and Kokoszka [70], and Hsing and Eubank [71], and a recent review is also available (Wang et al. [127]). The specific topic of amplitude and phase variation is discussed in [110, Chapter 7] and [127, Section 5.2]. The next paragraph gives some selective references.

One of the first functional registration techniques employed dynamic programming (Wang and Gasser [128]) and dates back to Sakoe and Chiba [118]. Landmark registration consists of identifying salient features for each curve, called landmarks, and aligning them (Gasser and Kneip [61]; Gervini and Gasser [63]). In pairwise synchronisation (Tang and Müller [122]) one aligns each pair of curves and then derives an estimator of the warp functions by linear averaging of the pairwise registration maps. Another class of methods involves a template curve, to which each observation is registered, minimising a discrepancy criterion; the template is then iteratively updated (Wang and Gasser [129]; Ramsay and Li [108]). James [72] defines a “feature function” for each curve and uses the moments of the feature function to guarantee identifiability. Elastic registration employs the Fisher–Rao metric that is invariant to warpings and calculates averages in the resulting quotient space (Tucker et al. [123]). Other techniques include semiparametric modelling (Rønn [115]; Gervini and Gasser [64]) and principal components registration (Kneip and Ramsay [82]). More details can be found in the review article by Marron et al. [90]. Wrobel et al. [131] have recently developed a registration method for functional data with a discrete flavour. It is also noteworthy that a version of the Wasserstein metric can also be used in the functional case (Chakraborty and Panaretos [34]).

The literature on the point processes case is more scarce; see the review by Wu and Srivastava [133].

A parametric version of Theorem 4.2.4 was first established by Bigot and Klein [22, Theorem 5.1] in , extended to a compact nonparametric formulation in Zemel and Panaretos [134]. There is an infinite-dimensional linear version in Masarotto et al. [91]. The current level of generality appears to be new.

Theorem 4.4.1 is a stronger version of Panaretos and Zemel [100, Theorem 1] where it was assumed that τ _n must diverge to infinity faster than . An analogous construction under the Bayesian paradigm can be found in Galasso et al. [58]. Optimality of the rates of convergence in Theorem 4.6.3 is discussed in detail by Bigot et al. [23], where finiteness of the functional J ₂ (see Sect. 4.7) is assumed and consequently is improved to .

As far as we know, Theorem 4.6.5 (taken from [100]) is the first central limit theorem for Fréchet means in Wasserstein space. When the measures Λ _i are observed exactly (no amplitude variation: τ _n = ∞ and σ = 0) Kroshnin et al. [84] have recently proven a central limit theorem for random Gaussian measures in arbitrary dimension, extending a previous result of Agueh and Carlier [3]. It seems likely that in a fully nonparametric setting, the rates of convergence (compare Theorem 4.6.3) might be slower than ; see Ahidar-Coutrix et al. [4].

The magnitude of the amplitude variation in Theorem 4.6.3 pertains to the rates of convergence of to zero (Sect. 4.7). This is a topic of intense research, dating back to the seminal paper by Dudley [46], where a version of Theorem 4.7.8 with p = 1 is shown for the bounded Lipschitz metric. The lower bounds proven in this section were adapted from [46], Fournier and Guillin [54], and Weed and Bach [130].

The version of Theorem 4.7.8 given here can be found in [130] and extends Boissard and Le Gouic [27]. Both papers [27, 130] work in a general setting of complete separable metric spaces. An additional term appears in the limiting case d = 2p, as already noted (for p = 1) by [46], and the classical work of Ajtai et al. [5] for μ uniform on [0, 1]². More general results are available in [54]. A longer (but far from being complete) bibliography is given in the recent review by Panaretos and Zemel [101, Subsection 3.3.1], including works by Barthe, Dobrić, Talagrand, and coauthors on almost sure results and deviation bounds for the empirical Wasserstein distance.

The J ₁ condition is due to del Barrio et al. [43], who showed it to be necessary and sufficient for the empirical process to converge in distribution to , with Brownian bridge. The extension to 1 ≤ p ≤∞ (and a lot more) can be found in Bobkov and Ledoux [25], employing order statistics and beta distributions to reduce to the uniform case. Alternatively, one may consult Mason [92], who uses weighted approximations to Brownian bridges.

An important aspect that was not covered here is that of statistical inference of the Wasserstein distance on the basis of the empirical measure. This is a challenging question and results by del Barrio, Munk, and coauthors are available for one-dimensional, elliptical, or discrete measures, as explained in [101, Section 3].

Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

References

1. 2.
   M. Agueh, G. Carlier, Barycenters in the Wasserstein space. Soc. Indus. Appl. Math. 43(2), 904–924 (2011)
2. 3.
   M. Agueh, G. Carlier, Vers un théorème de la limite centrale dans l’espace de Wasserstein? C.R. Math. 355(7), 812–818 (2017)
3. 4.
   A. Ahidar-Coutrix, T. Le Gouic, Q. Paris, Convergence rates for empirical barycenters in metric spaces: curvature, convexity and extendable geodesics. Probab. Theory Relat. Fields (2019). https://​link.​springer.​com/​article/​10.​1007%2Fs00440-019-00950-0
4. 5.
   M. Ajtai, J. Komlós, G. Tusnády, On optimal matchings. Combinatorica 4(4), 259–264 (1984)
5. 9.
   P.C. Álvarez-Esteban, E. del Barrio, J.A. Cuesta-Albertos, C. Matrán, A fixed-point approach to barycenters in Wasserstein space. J. Math. Anal. Appl. 441(2), 744–762 (2016)
6. 12.
   L. Ambrosio, N. Gigli, G. Savaré, Gradient Flows in Metric Spaces and in the Space of Probability Measures. Lectures in Mathematics. ETH Zürich, 2nd edn. (Springer, Berlin, 2008)
7. 14.
   E. Anderes, S. Borgwardt, J. Miller, Discrete Wasserstein barycenters: optimal transport for discrete data. Math. Meth. Oper. Res. 84(2), 1–21 (2016)
8. 16.
   A. Bandeira, P. Rigollet, J. Weed, Optimal rates of estimation for multi-reference alignment (2017). arXiv:1702.08546
9. 22.
   J. Bigot, T. Klein, Characterization of barycenters in the Wasserstein space by averaging optimal transport maps. ESAIM: Probab. Stat. 22, 35–57 (2018)
10. 23.
    J. Bigot, R. Gouet, T. Klein, A. López, Upper and lower risk bounds for estimating the Wasserstein barycenter of random measures on the real line. Electron. J. Stat. 12(2), 2253–2289 (2018)
11. 25.
    S. Bobkov, M. Ledoux, One-Dimensional Empirical Measures, Order Statistics and Kantorovich Transport Distances, vol. 261, no. 1259 (Memoirs of the American Mathematical Society, Providence, 2019). https://​doi.​org/​10.​1090/​memo/​1259
12. 27.
    E. Boissard, T. Le Gouic, On the mean speed of convergence of empirical and occupation measures in Wasserstein distance. Ann. Inst. H. Poincaré. Probab. Stat. 50(2), 539–563 (2014)
13. 34.
    A. Chakraborty, V.M. Panaretos, Functional registration and local variations: Identifiability, rank, and tuning (2017). arXiv:1702.03556
14. 41.
    D.J. Daley, D. Vere-Jones, An Introduction to the Theory of Point Processes: Volume II: General Theory and Structure (Springer, Berlin, 2007)
15. 43.
    E. del Barrio, E. Giné, C. Matrán, Central limit theorems for the Wasserstein distance between the empirical and the true distributions. Ann. Probab. 27(2), 1009–1071 (1999)
16. 46.
    R.M. Dudley, The speed of mean Glivenko–Cantelli convergence. Ann. Math. Stat. 40(1), 40–50 (1969)
17. 51.
    F. Ferraty, P. Vieu, Nonparametric Functional Data Analysis: Theory and Practice (Springer, Beriln, 2006)
18. 54.
    N. Fournier, A. Guillin, On the rate of convergence in Wasserstein distance of the empirical measure. Probab. Theory Rel. Fields 162(3–4), 707–738 (2015)
19. 58.
    B. Galasso, Y. Zemel, M. de Carvalho, Bayesian semiparametric modelling of phase-varying point processes (2018). arXiv:1812.09607
20. 61.
    T. Gasser, A. Kneip, Searching for structure in curve samples. J. Amer. Statist. Assoc. 90(432), 1179–1188 (1995)
21. 63.
    D. Gervini, T. Gasser, Self-modelling warping functions. J. Roy. Stat. Soc.: Ser. B 66(4), 959–971 (2004)
22. 64.
    D. Gervini, T. Gasser, Nonparametric maximum likelihood estimation of the structural mean of a sample of curves. Biometrika 92(4), 801–820 (2005)
23. 70.
    L. Horváth, P. Kokoszka, Inference for Functional Data with Applications, vol. 200 (Springer, Berlin, 2012)
24. 71.
    T. Hsing, R. Eubank, Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators (Wiley, Hoboken, 2015)
25. 72.
    G.M. James, Curve alignment by moments. Ann. Appl. Stat. 1(2), 480–501 (2007)
26. 73.
    H.E. Jones, N. Bayley, The Berkeley growth study. Child Dev. 12(2), 167–173 (1941)
27. 75.
    O. Kallenberg, Random Measures, 3rd edn. (Academic, Cambridge, 1983)
28. 79.
    A. Karr, Point Processes and Their Statistical Inference, vol. 7 (CRC Press, Boca Raton, 1991)
29. 82.
    A. Kneip, J.O. Ramsay, Combining registration and fitting for functional models. J. Amer. Stat. Assoc. 103(483), 1155–1165 (2008)
30. 84.
    A. Kroshnin, V. Spokoiny, A. Suvorikova, Statistical inference for Bures–Wasserstein barycenters (2019). arXiv:1901.00226
31. 88.
    E.L. Lehmann, A general concept of unbiasedness. Ann. Math. Stat. 22(4), 587–592 (1951)
32. 90.
    J.S. Marron, J.O. Ramsay, L.M. Sangalli, A. Srivastava, Functional data analysis of amplitude and phase variation. Stat. Sci. 30(4), 468–484 (2015)
33. 91.
    V. Masarotto, V.M. Panaretos, Y. Zemel, Procrustes metrics on covariance operators and optimal transportation of Gaussian processes. Sankhyā A 81, 172–213 (2019) (Invited Paper, Special Issue on Statistics on non-Euclidean Spaces and Manifolds)
34. 92.
    D.M. Mason, A weighted approximation approach to the study of the empirical Wasserstein distance, in ed. by C. Houdré, D.M. Mason, P. Reynaud-Bouret, J. Rosiński, High Dimensional Probability VII (Birkhäuser, Basel, 2016), pp. 137–154
35. 100.
    V.M. Panaretos, Y. Zemel, Amplitude and phase variation of point processes. Ann. Stat. 44(2), 771–812 (2016)
36. 101.
    V.M. Panaretos, Y. Zemel, Statistical aspects of Wasserstein distances. Annu. Rev. Stat. Appl. 6, 405–431 (2019)
37. 108.
    J. Ramsay, X. Li, Curve registration. J. Roy. Stat. Soc.: Ser. B (Stat. Methodol.) 60(2), 351–363 (1998)
38. 109.
    J.O. Ramsay, B.W. Silverman, Applied Functional Data Analysis: Methods and Case Studies, vol. 77 (Citeseer, 2002)
39. 110.
    J.O. Ramsay, B.W. Silverman, Functional Data Analysis, 2nd edn. (Springer, Berlin, 2005)
40. 111.
    J.O. Ramsay, H. Wickham, S. Graves, G. Hooker, FDA: Functional Data Analysis, R package version 2.4.8 (2018)
41. 115.
    B.B. Rønn, Nonparametric maximum likelihood estimation for shifted curves. J. Roy. Stat. Soc.: Ser. B (Stat. Methodol.) 63(2), 243–259 (2001)
42. 118.
    H. Sakoe, S. Chiba, Dynamic programming algorithm optimization for spoken word recognition. IEEE Trans. Acoust. Speech Signal Process. 26(1), 43–49 (1978)
43. 122.
    R. Tang, H.-G. Müller, Pairwise curve synchronization for functional data. Biometrika 95(4), 875–889 (2008)
44. 123.
    J.D. Tucker, W. Wu, A. Srivastava, Generative models for functional data using phase and amplitude separation. Comput. Stat. Data Anal. 61, 50–66 (2013)
45. 127.
    J.-L. Wang, J.-M. Chiou, H.-G. Müller, Functional data analysis. Annu. Rev. Stat. Appl. 3, 257–295 (2016)
46. 128.
    K. Wang, T. Gasser, Alignment of curves by dynamic time warping. Ann. Stat. 25(3), 1251–1276 (1997)
47. 129.
    K. Wang, T. Gasser, Synchronizing sample curves nonparametrically. Ann. Stat. 27, 439–460 (1999)
48. 130.
    J. Weed, F. Bach, Sharp asymptotic and finite-sample rates of convergence of empirical measures in Wasserstein distance. Bernoulli 25(4A), 2620–2648 (2019). https://​projecteuclid.​org/​euclid.​bj/​1568362038
49. 131.
    J. Wrobel, V. Zipunnikov, J. Schrack, J. Goldsmith, Registration for exponential family functional data. Biometrics 75, 48–57 (2019)
50. 132.
    S. Wu, H.-G. Müller, Z. Zhang, Functional data analysis for point processes with rare events. Stat. Sinica 23(1), 1–23 (2013)
51. 133.
    W. Wu, A. Srivastava, Analysis of spike train data: Alignment and comparisons using the extended Fisher–Rao metric. Electron. J. Stat. 8, 1776–1785 (2014)
52. 134.
    Y. Zemel, V.M. Panaretos, Fréchet means and Procrustes analysis in Wasserstein space. Bernoulli 25(2), 932–976 (2019)

Footnotes

1

For instance, the arithmetic average of two scalar Gaussians N(μ ₁, 1) and N(μ ₂, 1) will be their mixture with equal weights, but their Fréchet–Wasserstein average will be the Gaussian (see Lemma 4.2.1), which is arguably more representative from an intuitive point of view. In much the same way, the Fréchet–Wasserstein average of probability measures representing some type of object (e.g., normalised greyscale images of faces) will also be an object of the same type. This sort of phenomenon is well-known in manifold statistics, more generally, and is arguably one of the key motivations to account for the non-linear geometry of the sample space, rather than imbed it into a larger linear space and use the addition operation.

 

2

As the functional case will only serve as a motivation, our treatment of this case will mostly be heuristic and superficial. Rigorous proofs and more precise details can be found in the books by Ferraty and Vieu [51], Horváth and Kokoszka [70], or Hsing and Eubank [71]. The notion of amplitude and phase variation is discussed in the books by Ramsay and Silverman [109, 110] that are of a more applied flavour. One can also consult the review by Wang et al. [127], where amplitude and phase variation are discussed in Sect. 5.2.

 

3

If the cumulative count process Γ(t) = Π[0, t) is mean-square continuous, then the use of the term “amplitude variation” can be seen to remain natural, as Γ(t) will admit a Karhunen–Loève expansion, with all stochasticity being attributable to the random amplitudes in the expansion.

 

4

This can be defined as Bochner integral in the space of measurable bounded T : K → K.

 

5

In finite dimensional (or more generally, locally compact metric) spaces. If is an infinite-dimensional Hilbert space, the vague topology is trivial. This is stated and proved as Lemma 5 on page 27 in the supplement.