Wasserstein Covariance for Multiple Random Densities

TL;DR

This paper introduces a Wasserstein covariance measure for dependent density data, extending optimal transport-based analysis to multivariate and continuously-indexed densities, with theoretical guarantees and practical applications.

Contribution

It proposes a novel Wasserstein covariance framework for dependent densities, including estimators and asymptotic theory, applicable to neuroimaging and demographic data.

Findings

Wasserstein covariance effectively captures dependencies in density data.

Theoretical results establish estimator consistency and asymptotic normality.

Applications demonstrate utility in brain connectivity and mortality trends.

Abstract

A common feature of methods for analyzing samples of probability density functions is that they respect the geometry inherent to the space of densities. Once a metric is specified for this space, the Fr\'echet mean is typically used to quantify and visualize the average density from the sample. For one-dimensional densities, the Wasserstein metric is popular due to its theoretical appeal and interpretive value as an optimal transport metric, leading to the Wasserstein-Fr\'echet mean or barycenter as the mean density. We extend the existing methodology for samples of densities in two key directions. First, motivated by applications in neuroimaging, we consider dependent density data, where a $p$-vector of univariate random densities is observed for each sampling unit. Second, we introduce a Wasserstein covariance measure and propose intuitively appealing estimators for both fixed and diverging $p$, where the latter corresponds to continuously-indexed densities. We also give theory demonstrating consistency and asymptotic normality, while accounting for errors introduced in the unavoidable preparatory density estimation step. The utility of the Wasserstein covariance matrix is demonstrated through applications to functional connectivity in the brain using functional magnetic resonance imaging data and to the secular evolution of mortality for various countries.

Figures

Equations

$${\int }_{\mathbb{R}}{\left\{T^{\ast }(u)-u\right\}}^{2}f(u)\mathrm{d}u,$$

$$d_W^2(f,g)={\int }_{\mathbb{R}}{\left\{T(u)-u\right\}}^{2}f(u)\mathrm{d}u.$$

$$⟨T_1,T_2{⟩}_f={\int }_{\mathbb{R}}\left\{T_1(u)-u\right\}\left\{T_2(u)-u\right\}f(u)\mathrm{d}u,$$

$$f_{\oplus }=\underset{f\in \mathcal{D}}{\mathrm{argmin}}E\left\{d_W(\mathfrak{F},f{)}^2\right\},{\mathrm{var}}_{\oplus }(\mathfrak{F})=E\left\{d_W(\mathfrak{F},f_{\oplus }{)}^2\right\}.$$

$${\mathrm{var}}_{\oplus }({\mathfrak{F}}_j)=E\left\{d_W({\mathfrak{F}}_j,f_{\oplus ,j}{)}^2\right\}=E\left[{\int }_{\mathbb{R}}{\left\{T_j(u)-u\right\}}^2{f}_{\oplus ,j}(u)\mathrm{d}u\right]=E\left(⟨T_j,T_j{⟩}_{f_{\oplus ,j}}\right).$$

$${\mathrm{cov}}_{\oplus }({\mathfrak{F}}_1,{\mathfrak{F}}_2)=E\left(⟨{\tilde{T}}_1,T_2{⟩}_{f_{\oplus ,2}}\right)=E\left(⟨T_1,{\tilde{T}}_2{⟩}_{f_{\oplus ,1}}\right)={\mathrm{cov}}_{\oplus }({\mathfrak{F}}_2,{\mathfrak{F}}_1)$$

$$⟨T_1,T_2{⟩}_f={\int }_0^1\left\{F^{-1}(t)-G_1^{-1}(t)\right\}\left\{F^{-1}(t)-G_2^{-1}(t)\right\}\mathrm{d}t,$$

$${\mathrm{cov}}_{\oplus }({\mathfrak{F}}_j,{\mathfrak{F}}_k)=E\left[{\int }_0^1\left\{F_j^{-1}(t)-F_{\oplus ,j}^{-1}(t)\right\}\left\{F_k^{-1}(t)-F_{\oplus ,k}^{-1}(t)\right\}\mathrm{d}t\right](j,k=1,2).$$

$${\mathcal{C}}_{jk}(s,t)=\mathrm{cov}\left\{F_j^{-1}(s),F_k^{-1}(t)\right\}(j,k=1,2;0\le s,t\le 1),$$

$${\mathrm{cov}}_{\oplus }({\mathfrak{F}}_j,{\mathfrak{F}}_k)={\int }_0^1{\mathcal{C}}_{jk}(t,t)\mathrm{d}t=\text{Tr}({\mathcal{M}}_{jk}),$$

$${\left({\Sigma }_{\oplus }\right)}_{jk}={\mathrm{cov}}_{\oplus }({\mathfrak{F}}_j,{\mathfrak{F}}_k),$$

$${\Sigma }_{\oplus }(y,z)={\mathrm{cov}}_{\oplus }\left\{\mathcal{F}(y),\mathcal{F}(z)\right\}(0\le y,z\le 1).$$

$$\underset{f\in \mathcal{D}}{\sup }E\left\{d_W(f,\hat{f})\right\}=O(b_N).$$

$$\underset{f\in \mathcal{D}}{\sup }\underset{u\in D_f}{\sup }\max \left\{f(u),1/f(u),|f^{\prime }(u)|\right\}<\infty$$

$${\left({\hat{\Sigma }}_{\oplus }\right)}_{jk}={\int }_0^1{\hat{\mathcal{C}}}_{jk}(t,t)\mathrm{d}t.$$

$${\hat{\Sigma }}_{\oplus }(y_j,y_k)={\int }_0^1{\hat{C}}_{jk}(t,t)\mathrm{d}t,$$

$$|{\Sigma }_{\oplus }(y_1,z_1)-{\Sigma }_{\oplus }(y_2,z_2)|\le L_1(|y_1-y_2|+|z_1-z_2|)(0\le y_1,y_2,z_1,z_2\le 1).$$

$${\int }_{[0,1{]}^4}\mathrm{cov}\left\{X_1^c(s;y)X_1^c(s;z),X_1^c(t;y)X_1^c(t;z)\right\}\mathrm{d}s\mathrm{d}t\mathrm{d}y\mathrm{d}z<\infty .$$

$${\int }_0^1{\int }_0^1{\left\{{\Sigma }_{\oplus }(y,z)-{\hat{\Sigma }}_{\oplus }(y,z)\right\}}^2\mathrm{d}y\mathrm{d}z=O_p(n^{-1}+b_N^2).$$

$$S={\Vert \log {\hat{\Sigma }}_{\oplus }^n-\log {\hat{\Sigma }}_{\oplus }^m\Vert }_F^2+{\Vert \log {\hat{\Sigma }}_{\oplus }^n-\log {\hat{\Sigma }}_{\oplus }^a\Vert }_F^2+{\Vert \log {\hat{\Sigma }}_{\oplus }^m-\log {\hat{\Sigma }}_{\oplus }^a\Vert }_F^2.$$

Full text

1 Introduction

The analysis of samples of density functions or distributions is an important and challenging problem for modern statistical practice (Delicado, 2011). Examples include distributions of age at death or mortality for different countries, warping functions, or distributions of voxel-to-voxel correlations of functional magnetic resonance imaging signals for a sample of subjects. While functional principal component analysis using cross-sectional averaging can be directly applied for samples of density functions (Kneip & Utikal, 2001), more recently techniques have been developed that incorporate the geometric constraints inherent to the space of density functions. A popular metric for data where each data atom corresponds to a randomly sampled distribution or density is the Wasserstein metric, both for its theoretical appeal and its convincing empirical performance in various applications (Bolstad et al., 2003; Zhang & Müller, 2011; Bigot et al., 2016, 2017; Panaretos & Zemel, 2016).

We consider samples of density data where multiple random densities per subject are observed, that is, repeated realizations of a stochastic process defined on $\mathcal{D}^{p}$, where $\mathcal{D}$ is a space of one-dimensional density functions. For $p=1$, a variety of methods have been proposed, focusing on concepts of mean, modes of variation and dimensionality reduction (Delicado, 2011; Petersen & Müller, 2016; Hron et al., 2016). However, generalizations of these methods for $p>1$, where one has a vector of densities, have not yet been developed, even though such data arise in various applications. An essential first step is the extension of the mean and variance concepts for samples of densities to a vector of means and covariance matrix. We demonstrate here the usefulness of these concepts for applications. The type of data we consider is different from data where one observes, for each subject, a random sample of identically distributed $p$-vectors, for which the joint distribution is of interest. Instead, we focus on the joint modeling of $p$ univariate densities.

Using the Wasserstein geometry of optimal transport on the space of densities, we extend the concepts of Fréchet mean and variance to a measure of Wasserstein covariance between component densities for a $p$-variate density-valued process. As one does not observe the densities themselves but rather samples of univariate data that they generate, a preliminary density estimation step is necessary and is taken into account in our analysis. Our theoretical arguments show that the population Wasserstein covariance can be estimated consistently, with a limiting Gaussian distribution under sufficient conditions for the estimation error associated with the preliminary estimation step to be negligible. Motivated by one of the applications, we also consider continuously varying densities, where one observes a discretized version of a continuously evolving density-valued process, similar to densely observed repeated functional data (Park & Staicu, 2015; Chen et al., 2017). For this situation we target a Wasserstein covariance surface and develop theory for its consistent estimation while still accounting for errors in density estimation.

The utility of the proposed methodology is demonstrated through the analysis of functional magnetic resonance imaging and mortality data. In the application to brain imaging, we investigate differences in intra-regional functional connectivity between Alzheimer’s patients, cognitively normal subjects, and a third group of individuals diagnosed with mild cognitive impairment, a sign of increased risk for developing dementia. Here, intra-regional functional connectivity is quantified by the distribution of voxel-to-voxel correlations within a neighborhood, called a functional hub, corresponding to a small region in the brain. Multiple densities per subject are obtained for a number of hubs simultaneously, leading to a random vector of densities obtained for each subject in the sample. In the mortality application, we consider distributions of age at death over a range of calendar years, and compare the resulting Wasserstein covariance surfaces between a group of Eastern European countries and a second group of other nations.

2 Wasserstein Covariance

2.1 The Wasserstein Metric and Geometry

Let $\mathcal{D}$ be a class of one-dimensional densities such that $\int_{\mathbb{R}}u^{2}f(u){\rm d}u<\infty$ for all $f\in\mathcal{D}.$ For $f,\,g\in\mathcal{D}$, suppose $Y\sim f$ and consider the collection of nondecreasing maps $T^{*}:\mathbb{R}\rightarrow\mathbb{R}$, known as transports, such that $T^{*}(Y)\sim g.$ The optimal transport problem that leads to the Wasserstein metric is to find the transport that minimizes

[TABLE]

and the solution is known to be $T=G^{-1}\circ F$, where $F$ and $G$ are the distribution functions of $f$ and $g$, respectively (Ambrosio et al., 2008). The resulting squared Wasserstein distance is

[TABLE]

This metric arises from a local inner product (Ambrosio et al., 2008). If $g_{j}\in\mathcal{D}$ ($j=1,2)$ have distribution functions $G_{j},$ the optimal transports $T_{j}=G_{j}^{-1}\circ F$ for given $F$ reside in the tangent space of $f$, and for each $f$ an inner product between $T_{1}$ and $T_{2}$ can be defined by

[TABLE]

so that $d_{W}(f,g_{j})^{2}=\langle T_{j},T_{j}\rangle_{f}=\langle T_{j}^{-1},T_{j}^{-1}\rangle_{g_{j}}.$

2.2 Wasserstein Mean, Variance and Covariance

For a random density $\mathfrak{F}$ in $\mathcal{D}$, its Fréchet mean and Fréchet variance (Fréchet, 1948) are natural tools for relating the distributional properties of $\mathfrak{F}$ to the Wasserstein geometry of $\mathcal{D}$ in terms of first and second order behaviour. When the space $\mathcal{D}$ is endowed with the metric $d_{W}$, the Wasserstein–Fréchet, or simply Wasserstein, mean and variance of $\mathfrak{F}$ are

[TABLE]

For a single density process $\mathfrak{F},$ the theoretical and practical properties of the Wasserstein mean have been thoroughly investigated (Bolstad et al., 2003; Zhang & Müller, 2011; Panaretos & Zemel, 2016; Bigot et al., 2017), and recently the Wasserstein variance was adopted to quantify variability explained when performing dimension reduction for densities (Petersen & Müller, 2016). To quantify the dependence between two random densities, we propose here the extension of these concepts to a Wasserstein covariance measure.

For a bivariate density process $(\mathfrak{F}_{1},\mathfrak{F}_{2})$, where the $\mathfrak{F}_{j}$ are random elements of $\mathcal{D}$, $j=1,2,$ denote by $F_{\oplus,j}$ the distribution function of the Wasserstein mean of $\mathfrak{F}_{j}$ (4) and the corresponding density by $f_{\oplus,j}$. The random optimal transport from $f_{\oplus,j}$ to $\mathfrak{F}_{j}$ is $T_{j}=F_{j}^{-1}\circ F_{\oplus,j},$ where $F_{j}$ is the distribution function of $\mathfrak{F}_{j},$ so that the Wasserstein variances are

[TABLE]

This suggests defining a Wasserstein covariance measure as an expected inner product between $T_{1}$ and $T_{2}.$ Since in general $f_{\oplus,1}\neq f_{\oplus,2},$ these transports reside in different tangent spaces. Adopting a common technique for manifold-valued data (Yuan et al., 2012), we push $T_{1}$ to a new transport map $\tilde{T}_{1}$ in the tangent space of $f_{\oplus,2}$ by a parallel transport map. Intuitively, as $T_{1}$ is transported to $\tilde{T}_{1}$ along the geodesic connecting $f_{\oplus,1}$ to $f_{\oplus,2},$ its angle with the geodesic is preserved. In the Wasserstein geometry, setting $T_{\oplus,12}=F_{\oplus,1}^{-1}\circ F_{\oplus,2}$, one obtains $\tilde{T}_{1}=T_{1}\circ T_{\oplus,12}-T_{\oplus,12}+\textrm{id},$ where id is the identity map. Due to symmetry of this operation, if $\widetilde{T}_{2}=T_{2}\circ T_{\oplus,12}^{-1}-T_{\oplus,12}^{-1}+\textrm{id},$ then $\langle\widetilde{T}_{1},T_{2}\rangle_{f_{\oplus,2}}=\langle T_{1},\widetilde{T}_{2}\rangle_{f_{\oplus,1}}.$ Thus, the Wasserstein geometry motivates

[TABLE]

as the Wasserstein covariance between the two random densities.

2.3 Expression in Terms of Quantile Functions

The well known fact that for one-dimensional densities, the Wasserstein geometry is closely related to quantile functions (Villani, 2003; Panaretos & Zemel, 2016; Petersen & Müller, 2016) leads to a second characterization of Wasserstein covariance, which is useful for practice. The change of variable $t=F(u)$ applied to (3) gives the alternative expression

[TABLE]

so that $d_{W}^{2}(f,g)=\int_{0}^{1}\{F^{-1}(t)-G^{-1}(t)\}^{2}{\rm d}t.$ We can then express the quantities in (4) and (5) in terms of the random quantile functions $F_{j}^{-1}.$ The Wasserstein means $f_{\oplus,j}$ are characterized by their quantile functions $F_{\oplus,j}^{-1}(t)=E\{F_{j}^{-1}(t)\}$ ($0\leq t\leq 1;\,j=1,2),$ and the Wasserstein variances and covariances are

[TABLE]

Expressions (6), (7) reveal a connection to functional data analysis. Viewing $(F_{1}^{-1},F_{2}^{-1})$ as bivariate functional data, key objects are the mean functions $F_{\oplus,j}^{-1}(t)$ and covariance surfaces

[TABLE]

that characterize the first- and second-order behavior of the processes (Li & Hsing, 2010). Writing $\mathcal{M}_{jk}$ for the integral operator with kernel $\mathcal{C}_{jk},$ Fubini’s theorem implies that

[TABLE]

where $\textrm{Tr}(\cdot)$ is the operator trace. Accordingly, the Wasserstein variance of each component distribution can be interpreted as a summary of the variability in the quantile process, and Wasserstein covariance as a summary of covariability.

Using quantile functions has two major advantages when multiple densities are observed per subject. First, the derived notions of mean, variance, and covariance have geometric interpretations in the manifold induced by the Wasserstein metric in the space of distributions, as in (4) and (5). Second, quantile functions always have the same support $[0,1]$ regardless of the distributional supports of the $\mathfrak{F}_{j},$ so that the Wasserstein covariance remains well-defined even when the latter differ. In contrast, attempts to define similar covariance summaries based on cross-covariance operators of density or compositional representations are bound to fail when distributional supports differ. For example, if $\mathcal{G}_{jk}$ is the ordinary cross-covariance operator between densities $\mathfrak{F}_{j}$ and $\mathfrak{F}_{k}$, the operator trace is well-defined only when the supports coincide. Even when they are the same, taking $\textrm{Tr}(\mathcal{G}_{jk})$ as a summary covariance measure has no intuitive geometric meaning.

2.4 Wasserstein Covariance Matrices and Kernels

Considering a $p$-variate density process $\mathfrak{F}=(\mathfrak{F}_{1},\ldots,\mathfrak{F}_{p})$ with component Wasserstein means and variances $f_{\oplus,j}$ and ${\rm var}_{\oplus}(\mathfrak{F}_{j})$, $j=1,\ldots,p,$ we initially assume that $p$ is fixed, as in the brain imaging example in Section 4.1. The Wasserstein covariance matrix for $\mathfrak{F}$ is the $p\times p$ matrix with elements

[TABLE]

which is easily seen to be a valid covariance matrix.

Motivated by the mortality example in Section 4.2, suppose the components of $\mathfrak{F}$ are indexed by a continuous variable $0\leq y_{j}\leq 1$, which might represent time. To model the setting of repeatedly observed densities that are measured densely in time, we allow $p\rightarrow\infty$. Then $\mathfrak{F}$ can be thought of as a discretized version of an unobservable process $\mathcal{F}(y),$ $0\leq y\leq 1,$ with $\mathfrak{F}_{j}=\mathcal{F}(y_{j}).$ We target the Wasserstein mean surface $f_{\oplus}(\cdot;y)$ and covariance kernel

[TABLE]

3 Estimation of Wasserstein Covariance Objects

3.1 Density Estimation

While we defined the Wasserstein covariance for fully observed densities, in reality these densities are rarely if ever directly observed. Rather, the data actually available are collections of scalar random variables $W_{ijr}$ ($i=1,\ldots,n$; $j=1,\ldots p$; $r=1,\ldots,N_{ij}$), where $n$ is the number of subjects $i$, $p$ is the number of densities or distributions $j$ per subject, and $N_{ij}$ is the number of independent observations $r$ distributed according to the $j$-th density that are available for the $i$-th subject and may vary across $i,j$. The observed data can be viewed as resulting from two independent random mechanisms, where the first random mechanism generates the independent vectors of densities $f_{i}=\{f_{i1},\ldots,f_{ip}\}$ $(i=1,\ldots,n$), while the second generates the observations that are sampled from these distributions, $W_{ijr}\sim f_{ij}$. The $W_{ijr}$ are all independent and for each fixed $(i,j)$ are also identically distributed.

Obtaining density estimates $\hat{f}_{ij}$ from the observed data $W_{ijr}$ and using these as proxies for the $f_{ij}$, for the asymptotic analysis we need to address the challenge that these estimates are noisy and deviate from the true density targets. Since the targets $\Sigma_{\oplus}$ can be expressed as integrated moments of the multivariate quantile process, an obvious route would be to estimate the empirical quantile functions and proceed by averaging. However, this has some practical drawbacks and a preferred approach is to first construct a sample of smooth density estimates $\hat{f}_{ij}$, then obtaining smooth distribution functions by integration, quantile functions as inverse functions and the target quantities from the estimated quantile functions.

The theoretical analysis of the Wasserstein covariance estimates in Section 3.2 below requires the following assumption, where a preliminary density estimator is generically denoted by $\hat{f}$. {assumption}

There is a compact interval $I$ such that, for any $f\in\mathcal{D},$ its support $D_{f}$ is a compact interval contained in $I.$ If $W_{1},\ldots,W_{N}$ is a random sample from $f$, the density estimate $\hat{f}$ based on this sample is a probability density function on $D_{f}$ such that, for some decreasing sequence $b_{N}=o(1)$ as $N\rightarrow\infty,$

[TABLE]

A density estimator that satisfies Assumption 3.1 is described in Petersen & Müller (2016). If the support $D_{f}$ is known and the condition

[TABLE]

is satisfied, then one may take $b_{N}=N^{-1/3};$ see Proposition 1 in the Supplementary Material. With no assumed uniform lower bound on the random densities, Panaretos & Zemel (2016) proposed a density estimator for which $b_{N}=N^{-1/4}.$

Lastly, because $np$ densities need to be simultaneously estimated, each from data of varying sample sizes $N_{ij}$, these need to be tied to the number of independent subjects $n$, as follows. {assumption} There exists a sequence $N=N(n)$ such that $\min_{i,j}N_{ij}\geq N$ and $N\rightarrow\infty$ as $n\rightarrow\infty$.

Here $N$ is a uniform lower bound on the sample sizes for the $np$ densities to be estimated that must diverge with the number of subjects $n$ to ensure uniform consistency of the densities.

3.2 Wasserstein Covariance Estimation

For fixed $p$, given densities $f_{i}=(f_{i1},\ldots,f_{ip})$, $i=1,\ldots,n,$ that are independently and identically distributed according to $\mathfrak{F}=(\mathfrak{F}_{1},\ldots,\mathfrak{F}_{p})$, our main goal is to estimate the Wasserstein covariance matrix $\Sigma_{\oplus}$, with elements defined in (8). We compute density estimates $\hat{f}_{ij}$, which are then mapped to their quantile function estimates $\hat{X}_{ij}=\hat{F}_{ij}^{-1}$. Write $\hat{X}_{ij}^{c}=\hat{X}_{ij}-n^{-1}\sum_{i=1}^{n}\hat{X}_{ij}$ and $\hat{\mathcal{C}}_{jk}(s,t)=n^{-1}\sum_{i=1}^{n}\hat{X}_{ij}^{c}(s)\hat{X}_{ik}^{c}(t).$ To target the Wasserstein covariance $\Sigma_{\oplus},$ the results of Section 2.3 suggest the estimator

[TABLE]

Theorem 1 demonstrates the overall rate of convergence of the Wasserstein covariance estimator, establishing asymptotic normality when $N$ diverges sufficiently fast. While our focus is on the Wasserstein covariance, the same rate of convergence is also obtained for the full covariance surface $\hat{\mathcal{C}}_{jk}(s,t)$ as an estimator of $\mathcal{C}_{jk}(s,t)$; see Theorem 3 in the Supplementary Material, where also auxiliary results and proofs can be found.

Theorem 3.1**.**

Suppose Assumptions 3.1 and 3.1 hold, and that $\mathfrak{F}_{j}\in\mathcal{D}$ almost surely ($j=1,\ldots,p).$ Then $\lVert\Sigma_{\oplus}-\hat{\Sigma}_{\oplus}\rVert_{F}=O_{p}(n^{-1/2}+b_{N}),$ where $\lVert\cdot\rVert_{F}$ denotes the Frobenius norm. Additionally, if $n^{1/2}b_{N}$ converges to zero, then there exists a zero-mean $p\times p$ Gaussian matrix $\mathfrak{C}$ such that $n^{1/2}\left(\hat{\Sigma}_{\oplus}-\Sigma_{\oplus}\right)$ converges weakly to $\mathfrak{C}.$

The covariance structure of $\mathfrak{C}$ is a four-dimensional array defined in the Supplementary Material. Under regularity conditions, the density estimator in Petersen & Müller (2016) satisfies $b_{N}=N^{-1/3}$ so that asymptotic Gaussianity is obtained for $N=n^{q}$, $q>3/2.$ Faster rates of convergence $b_{N}=N^{-\rho}$, with $1/3<\rho<1/2$, can be obtained under additional smoothness conditions and for suitable density estimators, and then weaker conditions on $q$ will suffice.

In the continuously indexed case, the vectors of densites $f_{i}=(f_{i1},\ldots,f_{ip})$ correspond to discretized versions of independent and identically distributed realizations $\tilde{f}_{i}(\cdot\,;y)$ of a latent dynamic density surface $\mathcal{F}(y),$ $0\leq y\leq 1,$ where $f_{ij}$ are independently distributed for different $i$ as $\mathcal{F}(y_{ij}).$ For simplicity, we require {assumption}

The number of observation times $p=p(n)$ satisfies $np^{-1}=O(1)$, and these are equidistant and common for all subjects, i.e., $y_{ij}=(j-1)/(p-1),\,1\leq j\leq p.$

We then estimate $\Sigma_{\oplus}(y_{j},y_{k})$ by the sample Wasserstein covariance kernel estimator

[TABLE]

followed by linearly interpolating these discretized estimates to obtain $\hat{\Sigma}_{\oplus}(y,z)$ for any $0\leq y,z\leq 1.$ For this interpolation to be negligible, we require two additional assumptions. {assumption}

There exists a constant $L_{1}>0$ such that

[TABLE]

{assumption}

With $X_{i}(t\,;y)$ denoting the random quantile function corresponding to $\tilde{f}_{i}(\cdot\,;y)$ and $X_{i}^{c}(t\,;y)=X_{i}(t\,;y)-E\{X_{i}(t\,;y)\},$

[TABLE]

Theorem 3.2**.**

Suppose Assumptions 3.1–3.2 hold and that $\mathcal{F}(y)\in\mathcal{D}$ almost surely. Then

[TABLE]

4 Applications

4.1 Functional Connectivity in Brain Imaging

The study of functional connectivity in the brain involves the identification of voxels or regions which exhibit similar behaviour, as quantified by neuroimaging techniques such as electroencephalography and functional magnetic resonance imaging. Of special interest are connections that are present when subjects are in the resting state (Allen et al., 2014). Studies of connections between spatially remote regions revealed the so-called default mode network in the resting brain, including nodes of high centrality that have been characterized as functional connectivity hubs (Buckner et al., 2009). The strength of connections between neighboring voxels, as opposed to those between remote regions, contains important information related to various biological factors (Gao et al., 2016) and neurological diseases (Zalesky et al., 2012). The strength of local connections in a particular brain region has been quantified in various ways (Tomasi & Volkow, 2010; Zang et al., 2004).

Petersen & Müller (2016) demonstrated the utility of a relatively simple approach to quantifying local connectivity within functional connectivity hubs that uses probability density functions formed by smoothing histograms of pairwise temporal correlations between the signals of each voxel within a hub and the signal at its central seed voxel. We implement this approach for the $p=10$ hubs in Table 3 of Buckner et al. (2009). The resulting 10-dimensional vectors of densities that are obtained for each subject are then analyzed with the proposed Wasserstein covariance. The ten hubs of interest are located in the left and right parietal lobes, medial superior frontal lobe, medial prefrontal lobe, left and right middle frontal lobes, posterior cingulate/precuneus region, right supramarginal lobe and the left and right middle temporal lobes.

The densities were constructed for $171$ cognitively normal subjects, $65$ subjects with Alzheimer’s disease, and a third group of $86$ subjects with mild cognitive impairment, with details about data preprocessing as in Petersen et al. (2016). The densities were then used to compute estimated Wasserstein mean densities $\hat{f}_{\oplus,j}$ for each group separately, which are displayed for each of the 10 hubs in Figure 1. The various groups show remarkable similarities in terms of the average Wasserstein connectivity distributions across these ten functional hubs. In contrast, the Wasserstein covariance and correlation matrices in Figure 2 exhibit clear differences between the groups and provide valuable additional information. Correlations are overall stronger for the normal subjects and those with mild cognitive impairment and exhibit different patterns of dependency across hubs. Complete second order behaviour can be studied through the estimated quantile covariance function estimates $\hat{\mathcal{C}}_{jk}$, where for $p=10,$ there are 45 such functions for each group. We examine the off-diagonal elements for $j<k$ by plotting slices corresponding to the three quartiles $s=0.25,0.5,75$; see the Supplementary Material.

We can also test for group differences in Wasserstein covariance using bootstrap samples obtained by centering all quantile functions with respect to their group means. To construct a bootstrap sample under the null hypothesis that all groups have the same Wasserstein covariance matrix, we simply center each multivariate quantile process at its corresponding group Wasserstein mean. By pooling these centered processes, a bootstrap sample is obtained by sampling with replacement and then dividing into three groups of proper size. With estimates $\hat{\Sigma}_{\oplus}^{n}$, $\hat{\Sigma}_{\oplus}^{m}$ and $\hat{\Sigma}_{\oplus}^{a}$ for the groups of normal, mildly cognitively impaired, and Alzheimer’s subjects, the test statistic for the global null hypothesis that all groups are equal is

[TABLE]

Using 1000 bootstrap samples, the global $p$-value was found to be $p=0.015$. An alternative test, obtained by replacing the estimated Wasserstein covariance matrices with the corresponding correlations in the above statistic, resulted in $p=0.093.$

To further explore the differences between normal and Alzheimer’s groups, the matrices were reduced to a subset of rows/columns corresponding to the lateral hub pairs in the middle frontal, middle temporal and parietal regions, then reordered using Ward’s hierarchical clustering algorithm with three clusters, visualized in Figure 3. This demonstrates the presence of asymmetry in the Alzheimer’s disease group, which is absent in the normal group. Derflinger et al. (2011) and others report findings of similar asymmetries in the brains of Alzheimer’s patients.

One can ask how the Wasserstein covariance approach compares to established approaches in brain imaging. A commonly used measure of local connectivity is regional homogeneity (Zang et al., 2004), corresponding to Kendall’s coefficient of concordance between the BOLD signals at the seed voxel and those of its immediate neighbors. This scalar measure of regional homogenity can be computed for each hub and each subject, resulting in a sample of ordinary multivariate data. The group covariance and correlation matrices are shown in the Supplementary Material. Some common patterns are seen, most notably an increased number of negative correlations for the Alzheimer’s group. However, the differences are not as stark, and the regional homogeneity covariances do not reveal the asymmetry seen in the Wasserstein covariance analysis.

4.2 Distribution of Age at Death for Period Cohorts

To gain a better understanding of human longevity, the study of the temporal evolution of the distributions of age at death and their dependency structure over calendar time is of interest. The Human Mortality Database provides yearly mortality and population data for 38 countries at <www.mortality.org>, which have been previously analyzed with various functional data analysis techniques (Hyndman & Shang, 2010; Chiou & Müller, 2009).

For a given country and calendar year, the probability distribution for mortality can be represented by its density. Consider a country for which life tables are available for the years $y_{j}\,(j=1,\ldots,p)$. For integer-valued ages $a=0,\ldots,110$, the life table provides the size of the population $m_{a}$ at least $a$ years old, normalized so that $m_{0}=100000$. These life tables were converted to histograms of age-at-death, which we then smoothed, applying the hades package, available at http://www.stat.ucdavis.edu/hades/, with a smoothing bandwidth of $h=2$. This led to estimated densities of age-at-death on the age domain $[20\text{years},110\text{years}]$.

To illustrate the proposed methods for continuously varying densities, we considered 32 countries and identify a subgroup of $n_{E}=8$ countries located in Eastern Europe for comparison with the remaining $n_{O}=24$ countries. Densities were estimated for every year between 1985 and 2005. Wasserstein means are depicted in the Supplementary Material and are found to be quite similar for the two groups, demonstrating increasing longevity over calendar time. In contrast, the estimated Wasserstein covariance and correlation surfaces of these two groups differ quite drastically, as seen in Figure 4.

Examining the diagonals of the Wasserstein covariance plots, the Eastern European countries are characterized by stagnant Wasserstein variability until 1993 when it increases sharply, followed by a steady increase between 1998 and 2005. For the non-Eastern European countries, the Wasserstein variability is maximal in 1994. Wasserstein correlations reveal that mortality dependencies are high and roughly constant over time for the non-Eastern European countries, while they are weak for Eastern European countries between years before 1990 and those after 1995.

These exploratory findings are not altogether surprising, given the economic and societal upheaval in Eastern Europe beginning around 1990, which is seen to be reflected in the Wasserstein covariance.

In addition to exploratory analysis, we implemented a bootstrap test as described in the previous section. With estimated Wasserstein covariance surfaces $\hat{\Sigma}_{\oplus}^{E}$ and $\hat{\Sigma}_{\oplus}^{O}$, the test statistic was computed as the square root distance between the associated operators (Pigoli et al., 2014) and implemented by computing the Frobenius distance between the principal square roots of the discretized matrix estimates. Computing this statistic for 1000 bootstrap samples, the $p$-value for the difference between Wasserstein covariances of Eastern European and other countries was $0.287$, while it was $0.007$ using Wasserstein correlations, providing some evidence for differences in the distributions of the density processes.

5 Discussion

For studying the covariance structure of vectors of random densities, the Wasserstein approach is preferred over possible alternatives such as the transformation approach (Petersen & Müller, 2016) or compositional methods based on the Aitchison geometry (Egozcue et al., 2006). This is because of its convincing practical behavior for the construction of barycenters (Bolstad et al., 2003; Zhang & Müller, 2011) and the theoretically appealing connections with optimal transport.

Specifically, the transformation approach where densities are mapped to the entire Hilbert space $L^{2}$ by means of a suitable transformation such as the log-quantile density transformation, could be applied to the components of the vectors of densities, which would lead to unconstrained multivariate functional data. For such data, any one of numerous available measures of functional covariance and correlation (Leurgans et al., 1993; Dubin & Müller, 2005; Eubank & Hsing, 2008; Yang et al., 2011) could then be harnessed. One would need to choose among many possible covariance measures and transformation maps, none of which is isometric to the Wasserstein distance. The resulting metric distortions make such an approach difficult to interpret.

In contrast, the proposed Wasserstein covariance has a canonical interpretation as an expected value of inner products of optimal transport maps. This means that the proposed Wasserstein covariance is similar to the notion of a regular covariance for an appropriate inner product, and can be interpreted as a measure of the degree of synchronization of the movement of probability mass from the Fréchet means to the random components of a bivariate density process. It thus emerges as a natural and compelling extension of the Wasserstein–Fréchet variance. The quantile function representation of the Wasserstein covariance in (6) facilitates the joint Wasserstein analysis of $p$ one-dimensional distributions with Wasserstein covariance matrices and surfaces, enhancing the appeal of the proposed approach for practical applications.

While the geometric notion of covariance in (5) can be extended to the case where the sample densities have a multivariate domain, the implementation via quantile functions cannot be extended for this case, and therefore the asymptotic theory we provide remains limited to the case of one-dimensional densities. An extension to multivariate domains and analogous notions of covariance with respect to other metrics are open problems for future research.