Statistical inference for Bures-Wasserstein barycenters

TL;DR

This paper introduces the Bures-Wasserstein barycenter as a Fréchet mean of positive semi-definite operators, exploring its existence, uniqueness, and statistical properties with applications to quantum mechanics.

Contribution

It defines the Bures-Wasserstein barycenter for operators, establishes conditions for its existence and uniqueness, and analyzes its convergence and concentration properties in statistical inference.

Findings

Existence and uniqueness conditions for Bures-Wasserstein barycenters.

Convergence and concentration results for empirical barycenters.

Connections to optimal transportation and quantum mechanics applications.

Abstract

In this work we introduce the concept of Bures-Wasserstein barycenter $Q_*$, that is essentially a Fr\'echet mean of some distribution $\mathbb{P}$ supported on a subspace of positive semi-definite Hermitian operators $\mathbb{H}_{+}(d)$. We allow a barycenter to be restricted to some affine subspace of $\mathbb{H}_{+}(d)$ and provide conditions ensuring its existence and uniqueness. We also investigate convergence and concentration properties of an empirical counterpart of $Q_*$ in both Frobenius norm and Bures-Wasserstein distance, and explain, how obtained results are connected to optimal transportation theory and can be applied to statistical inference in quantum mechanics.

Equations

$$d_{BW}^2(Q,S)=\mathrm{tr}Q+\mathrm{tr}S-2\mathrm{tr}{\left(Q^{1/2}S{Q}^{1/2}\right)}^{1/2}.$$

$$\mathbb{P}\left(\mathcal{S}\cap {\mathbb{H}}_{++}(d)\right)>0,\mathbb{E}\mathrm{tr}S<+\infty .$$

$$\mathcal{V}(Q)\stackrel{\mathrm{def}}{=}{\int }_{\mathcal{S}}{d}_{BW}^2(Q,S)d\mathbb{P}(S).$$

$$Q_{\ast }\in \underset{Q\in {\mathbb{H}}_{+}(d)}{\mathrm{argmin}}\mathcal{V}(Q).$$

$$Q_{\ast }\in \underset{Q\in {\mathbb{H}}_{+}(d)\cap \mathbb{A}}{\mathrm{argmin}}\mathcal{V}(Q).$$

$${\mathcal{V}}_n(Q)\stackrel{\mathrm{def}}{=}\frac{1}{n}\sum _{i=1}^n{d}_{BW}^2(Q,S_i).$$

$$Q_n=\underset{Q\in {\mathbb{H}}_{+}(d)\cap \mathbb{A}}{\mathrm{argmin}}{\mathcal{V}}_n(Q),{\mathcal{V}}_n\stackrel{\mathrm{def}}{=}{\mathcal{V}}_n(Q_n).$$

$$\sqrt{n}\left(Q_n-Q_{\ast }\right)\rightharpoonup \mathcal{N}\left(0,\text{Upxi}\right),$$

$$\sqrt{n}{{\widehat{\text{Upxi}}}_{\mathbf{n}}}^{-1/2}\left(Q_n-Q_{\ast }\right)\rightharpoonup \mathcal{N}\left(0,\mathbf{Id}\right),$$

$$\mathcal{L}\left(\sqrt{n}{d}_{BW}(Q_n,Q_{\ast })\right)\rightharpoonup \mathcal{L}\left(\text{norm}{\xi }_F\right),$$

$$d_w\left(\mathcal{L}\left(\sqrt{n}{d}_{BW}(Q_n,Q_{\ast })\right),\mathcal{L}\left(\text{norm}{{\xi }_n}_F\right)\right)\to 0,$$

$$\text{norm}{Q_{\ast }^{-1/2}{Q}_n{Q}_{\ast }^{-1/2}-I}_F\le \frac{\mathtt{C}(d+t)}{\sqrt{n}},d_{BW}(Q_n,Q_{\ast })\le \frac{\mathtt{C}(d+t)}{\sqrt{n}},$$

$$\mathbb{P}\left\{\sqrt{\mathrm{tr}S}\ge t\right\}\le B{e}^{-b{t}^2}\text{for any}t\ge 0,$$

$$\sqrt{n}\left({\mathcal{V}}_n-{\mathcal{V}}_{\ast }\right)\rightharpoonup \mathcal{N}\left(0,\mathrm{Var}{d}_{BW}^2(Q_{\ast },S)\right),$$

$$\left|{\mathcal{V}}_n-{\mathcal{V}}_{\ast }\right|\le \max \left(\frac{c_1(B,b,t,d)}{\sqrt{n}},\frac{c_2(B,b,t,d)}{n}\right),$$

$$\rho ,\sigma \in {\mathbb{H}}_{+}(d),\mathrm{tr}\rho =1,\mathrm{tr}\sigma =1.$$

$$d_B^2(\rho ,\sigma )=2\left(1-{\mathcal{F}}^{1/2}(\rho ,\sigma )\right).$$

$$d_{W_2}^2\left(\mathcal{N}(0,Q),\mathcal{N}(0,S)\right)=\mathrm{tr}Q+\mathrm{tr}S-2\mathrm{tr}{\left(Q^{1/2}S{Q}^{1/2}\right)}^{1/2}.$$

$$d_{BW}^2(Q,S)=\text{norm}{\left(T_Q^S-I\right)Q^{1/2}}_F^2=\mathrm{tr}\left(T_Q^S-I\right)Q\left(T_Q^S-I\right),$$

$$T_Q^S$$

$$\displaystyle=S^{1/2}\bigl{(}S^{1/2}QS^{1/2}\bigr{)}^{-1/2}S^{1/2}=Q^{-1/2}\bigl{(}Q^{1/2}SQ^{1/2}\bigr{)}^{1/2}Q^{-1/2}.$$

$${\Pi }_{\mathbb{M}}\mathbb{E}{T}_Q^S={\Pi }_{\mathbb{M}}I,Q\in {\mathbb{H}}_{++}(d),$$

$$Q_{\ast }=\mathbb{E}{\left(Q_{\ast }^{1/2}S{Q}_{\ast }^{1/2}\right)}^{1/2}⟺\mathbb{E}{T}_{Q_{\ast }}^S=I.$$

$$\sqrt{n}\left(Q_n-Q_{\ast }\right)\rightharpoonup \mathcal{N}\left(0,\text{Upxi}\right)$$

$$\sqrt{n}{\widehat{\text{Upxi}}}_n^{-1/2}\left(Q_n-Q_{\ast }\right)\rightharpoonup \mathcal{N}\left(0,(\mathbf{Id}{)}_{\mathbb{M}}\right),$$

$$(\mathbf{A}{)}_{\mathbb{M}}:\mathbb{M}\to \mathbb{M},(\mathbf{A}{)}_{\mathbb{M}}(X)={\Pi }_{\mathbb{M}}\mathbf{A}(X),X\in \mathbb{M}.$$

$$T_{Q_n}^S=S^{1/2}{\left(S^{1/2}{Q}_n{S}^{1/2}\right)}^{-1/2}{S}^{1/2},T_{Q_n}^S\approx T_{Q_{\ast }}^S+{\mathbf{dT}}_{Q_{\ast }}^S(Q_n-Q_{\ast }),$$

$$\mathcal{L}\left(\sqrt{n}{d}_{BW}(Q_n,Q_{\ast })\right)\rightharpoonup \mathcal{L}\left(\text{norm}{Q_{\ast }^{1/2}{\mathbf{dT}}_{Q_{\ast }}^{Q_{\ast }}(Z)}_F\right),$$

$$d_{w}\Bigl{(}\mathcal{L}\left(\sqrt{n}d_{BW}(Q_{n},Q_{*})\right),\mathcal{L}\left(\norm{Q^{1/2}_{n}\bm{dT}_{Q_{n}}^{Q_{n}}(Z_{n})}_{F}\right)\Bigr{)}\rightarrow 0,$$

$$\mathcal{L}\left(\sqrt{n}{d}_{BW}(Q_n,Q_{\ast })\right)\rightharpoonup \mathcal{L}\left(\sum _{i,j=1}^d\frac{Z_{ij}^2}{2(q_i+q_j)}\right),$$

$$Q_n^{\prime }\stackrel{\mathrm{def}}{=}{Q}_{\ast }^{-1/2}{Q}_n{Q}_{\ast }^{-1/2}.$$

Full text

Bures-Wasserstein barycenter,

Wasserstein barycenter,

Hermitian operators,

Central limit theorem,

Concentration,

and

1 Introduction

Space of finite-dimensional Hermitian operators $\,\mathbb{H}(d)\,$ provides a powerful toolbox for data representation. For instance, in quantum mechanics it is used for mathematical description of physical properties of a quantum system, also known as observables. The reason is due to the fact that the measurements obtained in a physical experiment should be associated to real-valued quantities. Hermitian operators posses real-valued spectrum and satisfy the above requirement. A subspace $\,\mathrm{Sym}(d)\subset\mathbb{H}(d)\,$ of real-valued symmetric matrices is also of great interest: points in $\,\mathrm{Sym}(d)\,$ are widely used for description of systems in engineering applications, medical studies, neural sciences, evolutionary biology e.t.c. Usually such data sets are considered to be randomly sampled from an unknown distribution $\,\mathbb{P}\,$ (Goodnight and Schwartz (1997); Calsbeek and Goodnight (2009); Álvarez-Esteban et al. (2015); del Barrio et al. (2017); Gonzalez et al. (2017)), and statistical characteristics of $\,\mathbb{P}\,$ such as, in particular, mean and variance, appear to be of interest for further planning of an experiment and analysis of obtained results fur further development of natural science models. The current study focuses on a space of positive semi-definite Hermitian matrices $\,\mathbb{H}_{+}(d)\subset\mathbb{H}(d)\,$ and presents a possible approach to analysis and aggregation of relevant statistical information from data-sets, for which the linearity assumption might be violated. This makes classical Euclidean definitions of mean and variance not sensitive enough to capture effects of interest. This case appears extremely often in multiple contexts. As an example one can consider a data set that is represented as probability measures which belong to the same scale-location family, e.g. some astronomic measurements Alvarez-Esteban et al. (2018), Example 4.6. Non-linearity assumption requires the development of a novel toolbox suitable for further statistical analysis. In order to detect non-linear effects, we suggest to endow $\,\mathbb{H}_{+}(d)\,$ with the Bures-Wasserstein distance $\,d_{BW}\,$ which is recently introduced in a seminal paper Bhatia et al. (2018). It is defined as follows. For any pair of positive matrices $\,Q,S\in\mathbb{H}_{+}(d)\,$ it is written as:

[TABLE]

It is worth noting that being restricted to the space of symmetric positive definite matrices $\,\mathrm{Sym}_{++}(d)\,$, $\,d_{BW}\,$ boils down to a classical 2-Wasserstein distance between measures that belong to the same scale-location family (see e.g. Agueh and Carlier (2011), Section 6 or Álvarez-Esteban et al. (2015)). A more detailed discussion on this particular choice of the distance is presented in Section 2.1.

After choosing a proper distance, we are now ready to introduce a model of information aggregation the statistical properties of which are investigated in the current study. Let $\,\mathbb{P}\,$ be a probability distribution supported on some set $\,\mathcal{S}\subseteq\mathbb{H}_{+}(d)\,$. Further without loss of generality we assume, that $\,\mathbb{P}\,$ assigns positive probability to the intersection of $\,\mathcal{S}\,$ with space of positive definite Hermitian matrices $\,\mathbb{H}_{++}(d)\,$, and that the spectrum of its elemnts $\,S\in\mathcal{S}\,$ is on average bounded away from infinity:

Assumption 1**.**

[TABLE]

Two statistically important characteristics of $\,\mathbb{P}\,$ are Fréchet mean and Fréchet variance. The former one can be regarded as a typical representative of a data-set in hand, whereas the latter appears in analysis on data variability, see e.g. Del Barrio et al. (2015). We briefly recall both concepts below. For an arbitrary point $\,Q\in\mathbb{H}_{+}(d)\,$ Fréchet variance of $\,\mathbb{P}\,$ is defined as

[TABLE]

Classical Fréchet mean of $\,\mathbb{P}\,$ is a set of global minimizers of $\,\mathcal{V}(Q)\,$:

[TABLE]

However, in many cases we are interested in a minimizer, that belongs to some affine sub-space $\,\mathbb{A}\,$:

[TABLE]

For instance, such a necessity may arise while considering a random set of quantum density operators. For introduction to density operators theory one may look through Fano (1957). This example is considered in more details in Section 3.2. Note, that the setting (1.3) covers the setting (1.2). So without loss of generality we further address only (1.3). Obviously, the first crucial question concerns existence and uniqueness of $\,Q_{*}\,$. And positive answers on both issues, along with necessary conditions, are presented in Theorem 2.1. This immediately allows us to define the global Fréchet variance of $\,\mathbb{P}\,$ as $\,\mathcal{V}_{*}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}\mathcal{V}(Q_{*})\,$.

Given an i.i.d. sample $\,S_{1},...,S_{n}\,$, $\,S_{i}\overset{iid}{\sim}\mathbb{P}\,$, one constructs an empirical analogue of $\,\mathcal{V}(Q)\,$:

[TABLE]

An empirical Fréchet mean and global empirical variance also exist and unique:

[TABLE]

These facts follow from Theorem 2.1. This work studies convergence of the estimators $\,Q_{n}\,$ and $\,\mathcal{V}_{n}\,$ and investigate concentration properties of both quantities. The discussion of practical applicability of the obtained results is postponed to Section 3. There we explain their relation to optimal transportation theory and present a possible application to statistical analysis in quantum mechanics.

1.1 Contribution of the present study

Central limit theorem and

concentration of $\,Q_{n}\,$

The first main result of this study concerns asymptotic normality of the approximation error of population Fréchet mean by its empirical counterpart:

[TABLE]

where “$\,\rightharpoonup\,$” stands for weak convergence, and $\,\bm{\Upxi}\,$ is some covariance operator acting on the linear subspace $\,\mathbb{M}\subset\mathbb{H}(d)\,$ associated with affine subspace $\,\mathbb{A}\,$. From now on we use bold symbols e.g. $\,\bm{A},\bm{B}\,$, to denote operators, wheres classical ones i.e.$\,{A},{B}\,$ stand for matrices or vectors. This convergence result cannot be directly used for construction of asymptotic confidence sets because it relies on the unknown covariance matrix $\,\bm{\Upxi}\,$. However, Theorem 2.2 ensures, that this covariance matrix can be replaced by its empirical counterpart $\,\bm{\hat{\Upxi}_{n}}\,$:

[TABLE]

where $\,\bm{Id}\,$ denotes an identity operator. Along with asymptotic normality of $\,(Q_{n}-Q_{*})\,$, we are interested in the limiting distribution of $\,\mathcal{L}\left(\sqrt{n}d_{BW}(Q_{n},Q_{*})\right)\,$. Corollary 2.1 shows, that

[TABLE]

where $\,\xi\,$ is some normally distributed vector. Data-driven asymptotic confidence sets for $\,\sqrt{n}d_{BW}(Q_{n},Q_{*})\,$ are obtained by replacement of $\,\xi\,$ by its empirical counterpart $\,\xi_{n}\,$:

[TABLE]

where $\,d_{w}\,$ is a metric which induces weak convergence. Furthermore, we investigate concentration properties of $\,Q_{n}\,$ in both Frobenius norm and $\,d_{BW}\,$ metric. The following two bounds hold with h.p.:

[TABLE]

where $\,\mathtt{C}\,$ stands for some generic constant. A more detailed discussion is presented in Theorem 2.3 and Corollary 2.2 respectively. It is worth noting that concentration results are obtained under assumption of sub-Gaussianity of $\,S\,$ in the following sense:

Assumption 2** (Sub-Gaussianity of $\,\sqrt{\operatorname{tr}S}\,$).**

Let $\,\sqrt{\operatorname{tr}S}\,$ be sub-Gaussian, i.e.

[TABLE]

with some constants $\,B,b>0\,$.

All above-mentioned results are closely connected to convergence and concentration of empirical $\,2\,$-Wasserstein barycenters. For the sake of transparency we postpone this discussion to Section 3.1.

CLT and concentration for $\,\mathcal{V}_{n}\,$

We also show asymptotic normality of approximation error of $\,\mathcal{V}_{*}\,$ by $\,\mathcal{V}_{n}\,$ and prove concentration of $\,\mathcal{V}_{n}\,$:

[TABLE]

[TABLE]

where the latter result holds with h.p., and $\,c_{1}(B,b,t,d)\,$, $\,c_{2}(B,b,t,d)\,$ are constants which depend on sub-Gaussianity parameters $\,B,b\,$, dimension $\,d\,$, and parameter $\,t\,$. See Theorem 2.4 and Theorem 2.5 respectively.

The paper is organised as follows. Section 2 explains the obtained results in more details. Section 3 illustrates the connection to other scientific problems. Finally, Section 4 contains simulations and experiments on both artificial and real data-sets.

2 Results

This section presents obtained results in more details, and the first question we address is the particular choice of the distance.

2.1 Bures-Wasserstein distance

The original Bures metric appears in quantum mechanics in relation to fidelity measure between two quantum states and is used for measurement of quantum entanglement Marian and Marian (2008); Dajka et al. (2011). Let $\,\rho\,$, $\,\sigma\,$ be two quantum states. Mathematically speaking, this means that

[TABLE]

Fidelity of these states is defined as $\,\mathcal{F}(\rho,\sigma)=\left(\operatorname{tr}\sqrt{\rho^{1/2}\sigma\rho^{1/2}}\right)^{2}\,$. It quantifies “closeness” of $\,\rho\,$ and $\,\sigma\,$, see  Jozsa (1994). It is obvious, that in case of (2.1) Bures-Wasserstein distance turns into

[TABLE]

It is interesting to note, that the $\,d_{BW}\,$ distance appears not only one of the central distances, used in quantum mechanics, but also an object of extensive investigation in transportation theory Takatsu et al. (2011). Let $\,\mathscr{N}(0,Q)\,$ and $\,\mathscr{N}(0,S)\,$ be two centred Gaussian distributions. Then 2-Wasserstein distance between them is written as

[TABLE]

The case of Gaussian measures is naturally extended to measures that belong to a same scale-location family Alvarez-Esteban et al. (2018). In the last few years Wasserstein distance attracts a lot of attention of data scientists and machine learning community, as it takes into account geometrical similarities between objects, see e.g.  Gramfort et al. (2015); Flamary et al. (2018); Montavon et al. (2016). Due to this fact $\,d_{BW}\,$ satisfies the requirement of taking into account non-linearity of a data set under consideration. For more information on optimal transportation theory we recommend  Villani (2009).

Following Bhatia et al. (2018), we continue to investigate properties of $\,d_{BW}(Q,S)\,$. The next lemma presents an alternative analytical expression for the distance.

Lemma 2.1**.**

Let $\,Q,S\in\mathbb{H}_{+}(d)\,$ and $\,Q\succ 0\,$. Then (1.1) can be rewritten as

[TABLE]

where

[TABLE]

By \bigl{(}S^{1/2}QS^{1/2}\bigr{)}^{-1/2} we denote the pseudo-inverse matrix \left(\bigl{(}S^{1/2}QS^{1/2}\bigr{)}^{1/2}\right)^{+}.

Note, that in optimal transportation theory $\,T_{Q}^{S}\,$ is referred to as an optimal push-forward (optimal transportation map) between two centred normal distributions $\,\mathcal{N}(0,Q)\,$ and $\,\mathcal{N}(0,S)\,$. Following optimal transport notations it is denoted as $\,T^{S}_{Q}\#\mathcal{N}(0,Q)=\mathcal{N}(0,S)\,$.

For general notes on optimal transportation maps see Brenier (1991); for a particular case of scale-location and Gaussian families one may refer to Alvarez-Esteban et al. (2018); Takatsu et al. (2011). Lemma A.2 presents differentiability of the optimal map $\,T^{S}_{Q}\,$. It is one of the key-ingredients in the proof of main results of the present study. Note, that in case of $\,\mathbb{A}=\mathrm{Sym}_{++}(d)\,$ differentiability of $\,T^{S}_{Q}\,$ is obtained in Rippl et al. (2016). More technical details on properties of $\,d_{BW}\,$ are presented in Section A.2. However, for better understanding of the proofs of main results we highly recommend to at least look through Section A.1 which is dedicated to investigation of properties of $\,T^{S}_{Q}\,$ and its differential $\,\bm{dT}_{Q}^{S}\,$.

2.2 Existence and uniqueness of $\,Q_{*}\,$ and $\,Q_{n}\,$

Along with investigation of properties of the distance in hand and before moving to more general statistical questions, one should ask her- or himself, whether Fréchet mean $\,Q_{*}\,$ exists and, if so, is it unique or not? Let $\,\mathbb{M}\,$ be a linear subspace of $\,\mathbb{H}(d)\,$ associated to $\,\mathbb{A}\,$, i.e. the following representation holds: $\,\mathbb{A}=Q_{0}+\mathbb{M}\,$ for some $\,Q_{0}\in\mathbb{H}(d)\,$. We further assume that $\,\mathbb{A}\,$ has a non-empty intersection with the space of positive definite operators:

Assumption 3**.**

$\,\mathbb{H}_{++}(d)\cap\mathbb{A}\neq\emptyset\,$,

The next theorem ensures existence and uniqueness of the Fréchet mean (1.3).

Theorem 2.1** (Existence and uniqueness of Fréchet mean $\,Q_{*}\,$).**

Under Assumptions 1 and 3 there exists a unique positive-definite barycenter $\,Q_{*}\,$ of $\,\mathbb{P}\,$: $\,Q_{*}\succ 0\,$. Moreover, it is characterised as the unique solution of the equation

[TABLE]

where $\,\bm{\varPi_{\mathbb{M}}}\,$ is the orthogonal projector onto $\,\mathbb{M}\,$.

Note, that this result generalises the result for scale-location families in 2-Wasserstein space, presented in Álvarez-Esteban et al. (2015), Theorem 3.10 and originally obtained in a seminal work  Agueh and Carlier (2011), Theorem 6.1. Namely, if $\,\mathbb{A}=\mathrm{Sym}_{++}(d)\,$, then $\,Q_{*}\,$ exists, is unique, and is characterised as the unique solution of a fixed-point equation similar to (2.4)

[TABLE]

Existence, uniqueness, and measurability of the estimator $\,Q_{n}\,$ defined in (1.4) is a direct corollary of the above theorem. The proof of Theorem 2.1 is presented in Section A.3.

2.3 Convergence of $\,Q_{n}\,$ and $\,d_{BW}(Q_{n},Q_{*})\,$

Armed with the knowledge about properties of $\,d_{BW}\,$, $\,Q_{*}\,$, and $\,Q_{n}\,$, we are now equipped enough, so that to introduce the main results of the current study. Theorem 2.2 presents asymptotic convergence of $\,Q_{n}\,$ to $\,Q_{*}\,$.

Theorem 2.2** (Central limit theorem for the Fréchet mean).**

Under Assumptions 1 and 3 an approximation error rate of the Fréchet mean $\,Q_{*}\,$ by its empirical counterpart $\,Q_{n}\,$ is

[TABLE]

where $\,\bm{\Upxi}\,$ is a self-adjoint linear operator acting from $\,\mathbb{M}\,$ to $\,\mathbb{M}\,$ defined in (A.7). Moreover, if $\,\operatorname{Var}\left(T^{S}_{Q_{*}}\right)\,$ is non-degenerated, then

[TABLE]

with $\,\bm{\hat{\Upxi}_{n}}\,$ is a data-driven empirical counterpart of $\,\bm{\Upxi}\,$ defined in (A.8).

Remark 1**.**

Here $\,(\bm{A})_{\mathbb{M}}\,$ denotes a restriction of a quadratic form $\,\bm{A}\,$ to a subspace $\,\mathbb{M}\,$:

[TABLE]

We intentionally postpone the explicit definitions of $\,\bm{\Upxi}\,$ and $\,\bm{\hat{\Upxi}_{n}}\,$, as they require an introduction of many technical details. This would make the description of main results less transparent. The proof of the theorem relies on the Fréchet differentiablilty of $\,T_{Q}^{S}\,$ in the vicinity of $\,Q_{*}\,$:

[TABLE]

where $\,\bm{dT}_{Q_{*}}^{S}\,$ is a differential of $\,T_{Q}^{S}\,$ at point $\,Q_{*}\,$. Here we imply differentiability of $\,T_{Q}^{S}\,$ by the lower argument $\,Q\,$.

It is worth noting that the result (B) obtained in CLT enables construction of data-driven asymptotic confidence sets. However, there might appear technical problems with inversion of the empirical covariance. For instance, numerical simulations show, that $\,\bm{\Upxi}\,$ can be degenerated if $\,\mathbb{P}\,$ is supported on a set of diagonal matrices. This immediately raises a question concerning the development of some other confidence set construction methodology based on re-sampling techniques which would simplify the process from computational point of view. We consider this as a subject for further research.

As soon as the Bures-Wasserstein distance is the main tool for the analysis in $\,\mathbb{H}_{+}(d)\,$, the convergence properties of $\,d_{BW}(Q_{n},Q_{*})\,$ are also of great interest. The next lemma is almost a straightforward corollary of Theorem 2.2.

Corollary 2.1** (Asymptotic distribution of $\,d_{BW}(Q_{n},Q_{*})\,$).**

Under conditions of Theorem 2.2 it holds

[TABLE]

where $\,Z\in\mathbb{M}\subset\mathbb{H}(d)\,$ is random matrix, s.t. $\,Z\sim\mathscr{N}\left(0,\bm{\Upxi}\right)\,$. Moreover,

[TABLE]

where $\,Z_{n}\in\mathbb{M}\,$ and $\,Z_{n}\sim\mathscr{N}\left(0,\bm{\Upxi_{n}}\right)\,$.

To illustrate the result, we consider the case of diagonal $\,Q_{*}\,$. This setting allows us to write down the explicit form of the limiting distribution. If $\,Q_{*}=\operatorname{diag}(q_{1},...,q_{d})\,$, then right-hand side of the above corollary for $\,Z\,$-case is:

[TABLE]

where $\,(Z_{ij})^{d}_{i,j=1}=Z\,$. All proofs are collected in Section A.3. Section 4 illustrates asymptotic behaviour of $\,\mathcal{L}\left(\sqrt{n}\norm{Q_{n}-Q_{*}}_{F}\right)\,$ and $\,\mathcal{L}\left(\sqrt{n}d_{BW}(Q_{n},Q_{*})\right)\,$ on both artificial and real data sets.

2.4 Concentration of $\,Q_{n}\,$

The next important issue is concentration properties of $\,Q_{n}\,$ under the assumption of sub-Gaussianity of $\,\mathbb{P}\,$ (Assumption 2).

Theorem 2.3** (Concentration of $\,Q_{n}^{\prime}\,$).**

Let

[TABLE]

It holds under Assumptions 2 and 3, that

[TABLE]

for any $\,t\geq 0\,$ and $\,n\geq n_{0}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}c_{Q}^{2}(d+t)^{2}\,$, where

[TABLE]

operator $\,\bm{F^{\prime}}\,$ is defined in (B.3), $\,\sigma_{T}\,$ comes from Proposition 1, and $\,\sigma_{F}\,$ and $\,U\,$ are defined in Lemma B.3.

Concentration of $\,d_{BW}(Q_{n},Q_{*})\,$ is a corollary of the above theorem.

Corollary 2.2** (Concentration of $\,Q_{n}\,$ in $\,d_{BW}\,$ distance).**

Under conditions of Theorem 2.3 the following result holds

[TABLE]

Proofs are collected in Section B.

2.5 Central limit theorem and concentration for $\,\mathcal{V}_{n}\,$

In this section we investigate properties of the Fréchet variance $\,\mathcal{V}_{n}\,$, defined in (1.4). The next theorem presents central limit theorem for empirical variance $\,\mathcal{V}_{n}\,$.

Theorem 2.4** (Central limit theorem for $\,\mathcal{V}_{n}\,$).**

Let $\,\mathbb{P}\,$ be s.t. $\,\mathbb{E}(\operatorname{tr}S)^{2}<\infty\,$ and $\,\mathbb{P}\left(\mathbb{H}_{++}(d)\right)=1\,$. Then

[TABLE]

The last important result of the current study describes concentration properties of $\,\mathcal{V}_{n}\,$.

Theorem 2.5** (Concentration of $\,\mathcal{V}_{n}\,$).**

Let Assumption 2 be fulfilled. Then under conditions of Theorem 2.4 it holds:

[TABLE]

with

[TABLE]

There$\,(\nu,b)\,$ are parameters of sub-exponential r.v. $\,d^{2}_{BW}(Q_{*},S)\,$.

Proofs of these two theorems are collected in Section B.1.

3 Connection to other problems

In this section we explain the connection of obtained results to some other problems. Section 3.1 investigates the relation between Bures-Wasserstein barycenter and 2-Wasserstein barycenter of some scale-location family. Section 3.2 illustrates the idea of search of a barycenter on some affine subspace $\,\mathbb{A}\subset\mathbb{H}(d)\,$.

3.1 Connection to scale-location families of measures

We first present the concept of a scale-location family of absolutely continuous measures supported on $\,\mathbb{R}^{d}\,$.

Definition 3.1**.**

Let $\,X\sim\mu\,$ be a random variable that follows law $\,\mu\,$: $\,\mu\in\mathcal{P}^{ac}_{2}(\mathbb{R}^{d})\,$, where $\,\mathcal{P}^{ac}_{2}(\mathbb{R}^{d})\,$ is a set of all continuous measures with finite second moment. A set of all affine transformations of $\,X\,$ is

[TABLE]

It is referred to as a scale-location family.

Scale-location families attract lots of attention in modern data analysis and appear in many practical applications, as this concept is user-friendly in terms of theoretical analysis and, at the same time, possess very high modelling power. For example, it is widely used in medical imaging Wassermann et al. (2010), modelling of molecular dynamic Gonzalez et al. (2017), clustering procedures del Barrio et al. (2017), climate modelling Mallasto and Feragen (2017), embedding of complex objects in low dimensional spaces Muzellec and Cuturi (2018) and so on.

A possible metric that takes into account non-linearity of the underlying data-set is $\,2\,$-Wasserstein distance $\,d_{W_{2}}\,$. Let $\,\mu_{1},\mu_{2}\,$ be elements of $\,\mathcal{SL}(\mu)\,$ and let $\,X\sim\mu_{1}\,$, $\,Y\sim\mu_{2}\,$. We denote their first and second moments as

[TABLE]

It is a well-known fact, that in case of scale-location families $\,d_{W_{2}}\,$ depends only on the first and second moments of observed measures:

[TABLE]

For more details on general class of optimal transportation distances we recommend excellent books Ambrosio and Gigli (2013) or Villani (2009).

Distribution over scale-location family

In many cases we are interested in scale-location families generated at random. Let $\,\mathbb{P}\,$ be a probability measure supported on some $\,\mathcal{SL}(\mu)\,$. And let $\,\left(\varOmega,F,\mathbb{P}\right)\,$ be a generic probability space, s.t. for any $\,\omega\in\varOmega\,$ there exists an image in \,\mu_{\omega}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}\mathcal{L}\bigl{(}P_{\omega}X+p_{\omega}\bigr{)}\,, where $\,P_{\omega}\in\mathrm{Sym}_{+}(d)\,$ is a scaling parameter and $\,p_{\omega}\in\mathbb{R}^{d}\,$ is a shift parameter. A randomly sampled measure $\,\mu_{\omega}\,$ belongs to $\,\mathcal{SL}(\mu)\,$ by construction, and its first and second moments $\,(m_{\omega},S_{\omega})\,$ are written as

[TABLE]

where $\,(r,Q)\,$ denote the first and the second moments of $\,\mu\,$.

Fréchet variance of $\,\mathbb{P}\,$ at any arbitrary point $\,\mu\,$ is written as

[TABLE]

Given an i.i.d. sample $\,\nu_{1},...,\nu_{n}\,$ from $\,\mathbb{P}\,$, we define an empirical analogon of $\,\mathcal{V}\,$:

[TABLE]

Then population and empirical barycenters $\,\mu_{*}\,$ and $\,\mu_{n}\,$ are

[TABLE]

Note, that $\,\mu_{*}\,$ and $\,\mu_{n}\,$ belong to $\,\mathcal{SL}(\mu)\,$ and are uniquely characterised by their first and second moments $\,(r_{*},Q_{*})\,$ and $\,(r_{n},Q_{n})\,$ respectively, see e.g. Theorem 3.10 Álvarez-Esteban et al. (2015):

[TABLE]

[TABLE]

It is worth noting that the concept of Wasserstein barycenter originally presented in a seminal work by  Agueh and Carlier (2011) becomes a topic of extensive scientific interest in the last few years. A work Bigot et al. (2012) focuses on convergence of parametric class of barycenters, while Bigot et al. (2017) investigate asymptotic properties of regularised barycenters. The most general results on limiting distribution of convergence of empirical barycenters are obtained in Ahidar-Coutrix et al. (2018). This work provides rates of convergence for empirical barycentres of a Borel probability measure on a metric space either under assumptions on weak curvature constraint of the underlying space or for a case of a non-negatively curved space on which geodesics, emanating from a barycenter, can be extended. Theorem 2.2 specifies the results, obtained in Ahidar-Coutrix et al. (2018) for the case of scale-location families. Corollary 2.1 partially answers an (implicit) question, raised by work

Le Gouic and Loubes (2015), concerned the rate of convergence of $\,d_{W_{2}}(\mu_{n},\mu_{*})\,$. Namely, for the case of scale-location families it is of order $\,\frac{1}{\sqrt{n}}\,$. However, the above mentioned work covers only (1.2) case. The paper Kroshnin (2018) obtains an analog of law of large numbers for the case of arbitrary cost functions for barycenters on some affine sub-space $\,\mathbb{A}\,$ (1.3). A result, similar in spirit to Theorem 2.5 is obtained in Del Barrio et al. (2016). However, there authors consider only the space of probability measures supported on the real line (i.e. $\,d=1\,$) endowed with 2-Wasserstein distance. To the best of our knowledge, there are no results similar to concentration Theorem 2.3 and Theorem 2.2 in case of 2-Wasserstein distance.

3.2 Connection to quantum mechanics

This section illustrates the idea of barycenter restricted to some affine sub-space $\,\mathbb{A}\,$. We first briefly recall the concept of quantum densities. Quantum density operator is used in quantum mechanics as a possible way of description of statistical state of a quantum system. It might be considered as an analogue to a phase-space density in classical statistical mechanics. The formalism was introduced by John von Neumann in 1927. In essence a density matrix $\,\rho\,$ is a Hermitian positive semi-definite operator with the unit trace, $\,\rho\in\mathbb{H}_{+}(d)\,$, $\,\operatorname{tr}\rho=1\,$.

Given a random ensemble of density matrices, one is able to recovery the mean using averaging in classical Euclidean sense. However, Bures-Wasserstein barycenter suggests an alternative way to define the “most typical” representant (1.3) in terms of fidelity measure (2.2). We consider a following statistical setting. Let $\,(\varOmega,F,\mathbb{P})\,$ be some mechanism which generates quantum states $\,\rho_{\omega}\,$. Given an i.i.d sample $\,\rho_{1},...,\rho_{n}\,$ we write population and empirical variance of $\,\mathbb{P}\,$ as

[TABLE]

Then population and empirical barycenters in the class of all $\,d\times d\,$-dimensional density operators are defined as

[TABLE]

It can be easily shown, that “taking global Fréchet barycenter” or, in other words neglecting the condition $\,\sigma:~{}\operatorname{tr}\sigma=1\,$, we end up with the global baryceneter, which is the solution of the fixed point equation which is already mentioned in Section 2: $\,\rho=\int\left(\rho^{1/2}\rho_{\omega}\rho^{1/2}\right)^{1/2}\mathbb{P}(d\omega)\,$. This is a contraction mapping. Thus $\,\operatorname{tr}\rho_{*}<1\,$ and $\,\rho_{*}\,$ is not a density operator. In other words condition $\,\operatorname{tr}\sigma=1\,$ ensures, that $\,\rho_{*}\,$ and $\,\rho_{n}\,$ also belong to the class of density operators. Taking into account the results obtained in Section 2, $\,\rho_{n}\,$ is a natural consistent estimator of $\,\rho_{*}\,$ with known rate of convergence and deviation properties.

4 Experiments on simulated and real data sets

4.1 Simulated data

In this section we consider a simulated data set. So as to generate a covariance matrix $\,S=U^{*}\varLambda U\,$, $\,S\in\mathrm{Sym}_{++}(d)\,$ we generate at random an orthogonal matrix $\,U\,$ and a diagonal matrix $\,\varLambda=\text{diag}(\lambda_{1},...,\lambda_{d})\,$, s.t. $\,\lambda_{i}\sim\text{Unif}[18,22]\,$. The following images Fig.1 - Fig.3 illustrate convergence of $\,\mathcal{L}\left(\sqrt{n}\|Q_{n}-Q_{*}\|_{F}\right)\,$, $\,\mathcal{L}\left(\sqrt{n}d_{BW}(Q_{n},Q_{*})\right)\,$, and $\,\mathcal{L}\left(\sqrt{n}(\mathcal{V}_{*}-\mathcal{V}_{n})\right)\,$ presented in Theorem 2.2, Corollary 2.1, and Theorem 2.4 respectively. The following numerical experiments were performed using  R . The population barycenter $\,Q_{*}\,$ was computed using a sample of $\,20000\,$ observed covariance matrices. A solid line depicts the density of a limiting distribution, whereas dashed lines correspond to densities for different sample sizes for Bures-Wasserstein barycenter $\,Q_{n}\,$ with $\,n=\{3,10,100,1000\}\,$. Simulation were carried out for matrices of size $\,d=5\,$ and $\,d=10\,$.

4.2 Data aggregation in climate modelling

In this section we carry out the experiments on a family of Gaussian process, using a climate-related data set, collected in Siberia (Russia) between 1930 and 2009 Bulygina and Razuvaev (2012); Tatusko (1990). We set $\,\mathcal{SL}(\mu)\,$ to be a family of Gaussian curves, that describe the daily minimum temperatures within one year, measured at a set of 30 randomly sampled meteorological stations. Each curve is obtained by means of regression and maximum likelihood estimation and is sampled in $\,50\,$ points. More details on this data set are provided in Mallasto and Feragen (2017). The scale-location family under consideration is written as

[TABLE]

where $\,\nu_{t}=GP(m_{t},S_{t})\,$ is a Gaussian process, characterised by mean $\,m_{t}\,$ and covariance $\,S_{t}\,$ inherent to a year $\,t\,$, $\,t\in\{1933,...,2009\}\,$. We let $\,m_{t}=0\,$ for all $\,t\,$. A Gaussian process $\,\mu_{*}=GP(r_{*},Q_{*})\,$ is the population Wasserstein barycenter of $\,\mathcal{SL}(\mu)\,$. It is characterised by $\,(r_{*},Q_{*})\,$

[TABLE]

A family of approximating processes$\,\mu_{n}=GP(r_{n},Q_{n})\,$ with parameters $\,(r_{n},Q_{n})\,$ (3.3) is constructed by means of re-sampling with replacement of the original data set. Sample size $\,n\,$ varies in range $\,n=\{3,10,70\}\,$. Fig. 5 and Fig. 5 present densities of $\,\mathcal{L}\left(\sqrt{n}\norm{Q_{n}-Q_{*}}_{F}\right)\,$ and $\,\mathcal{L}\left(\sqrt{n}d_{BW}(Q_{n},Q_{*})\right)\,$ respectively.

Appendix A Proof of Central Limit Theorem

A.1 Properties of $T_{Q}^{S}$

Proof of Lemma 2.1.

First, we prove that optimal $\,T\,$ is self-adjoint. Indeed, assume the opposite, then

[TABLE]

and thus $\,\operatorname{tr}Q^{1/2}TQ^{1/2}<\operatorname{tr}\left(Q^{1/2}SQ^{1/2}\right)^{1/2}\,$. Therefore

[TABLE]

If $\,T\,$ is Hermitian but not positive semi-definite, then $\,Q^{1/2}TQ^{1/2}\preccurlyeq\left(Q^{1/2}SQ^{1/2}\right)^{1/2}\,$, $\,Q^{1/2}TQ^{1/2}\neq\left(Q^{1/2}SQ^{1/2}\right)^{1/2}\,$, hence again $\,\operatorname{tr}Q^{1/2}TQ^{1/2}<\operatorname{tr}\left(Q^{1/2}SQ^{1/2}\right)^{1/2}\,$.

Finally, if $\,T\in\mathbb{H}_{+}(n)\,$, then it is straightforward to check that $\,T=T_{Q}^{S}\,$ given by (2.3) and

[TABLE]

∎

The proof of the Central Limit theorem mainly relies on the differentiability of the map (2.3). Lemma A.2 shows that $\,T_{Q}^{S}\,$ can be linearised in the vicinity of $\,Q\,$:

[TABLE]

where $\,\bm{dT}_{Q}^{S}\colon\mathbb{H}(d)\to\mathbb{H}(d)\,$ is a self-adjoint negative-definite operator and $\,\norm{\cdot}\,$ stands for an operator norm. Properties of $\,\bm{dT}_{Q}^{S}\,$ are investigated in Lemma A.3. Let us introduce some notation: if $\,G(A)\,$ is a functional of a matrix $\,A\,$, then we denote its differential as $\,\bm{d_{A}}G\,$.

Lemma A.1**.**

Map $\,Q\mapsto g(Q)=Q^{1/2}\,$ is differentiable on $\,\mathbb{H}_{++}(d)\,$, and its differential is given by

[TABLE]

where $\,Q=U^{*}\operatorname{diag}(q)U\,$ is the eigenvalue decomposition.

Proof.

First, let us consider the map $\,P\mapsto f(P)=P^{2}\,$. It is smooth and its differential

[TABLE]

is non-degenerated:

[TABLE]

whenever $\,P\in\mathbb{H}_{++}(d)\,$. From now on $\,\langle\cdot,\cdot\rangle\,$ denotes a scalar product associated to Frobenius norm.

Now applying the inverse function theorem we obtain that the inverse map $\,g(\cdot)\,$ is also smooth and its differential enjoys the following equation

[TABLE]

thus

[TABLE]

[TABLE]

and

[TABLE]

∎

Lemma A.2** (Fréchet-differentiability of the map $\,T^{S}_{Q}\,$).**

For any $\,S\in\mathbb{H}_{+}(d)\,$ the map $\,T_{Q}^{S}\,$ can be linearised in the vicinity of $\,Q\in\mathbb{H}_{++}(d)\,$ as

[TABLE]

where

[TABLE]

$\,U^{*}\varLambda U\,$* is an eigenvalue decomposition of $\,S^{1/2}QS^{1/2}\,$*

[TABLE]

Proof.

The proof mainly relies on the differentiation of the pseudo-inverse term \,\bigl{(}S^{1/2}QS^{1/2}\bigr{)}^{-1/2}\,, as soon as

[TABLE]

Obviously we can consider only restriction to $\operatorname{range}(S)$ and therefore assume w.l.o.g. $S\succ 0$. As $\,\left(S^{1/2}(Q+X)S^{1/2}\right)^{-1/2}=U^{*}\left(\varLambda+\varDelta\right)^{-1/2}U\,$, by Lemma A.1 and von Neumann series expansion we obtain for infinitesimal $\,X\in\mathbb{H}(d)\,$ and corresponding $\,\varDelta\,$ that

[TABLE]

Then the differential $\,\bm{d_{Q}}\left(S^{1/2}QS^{1/2}\right)^{-1/2}(X)\,$ is written as

[TABLE]

Therefore,

[TABLE]

where $\,\bm{dT}_{Q}^{S}(X)\,$ is defined by (A.1). ∎

Lemmas A.3 and A.4 are technical and explore properties of $\,\bm{dT}_{Q}^{S}\,$.

Lemma A.3**.**

For any $\,S\in\mathbb{H}_{+}(d)\,$, $\,Q\in\mathbb{H}_{++}(d)\,$, the properties of operator $\,\bm{dT}_{Q}^{S}\,$ defined in (A.1) are following:

1. (I)

it is self-adjoint; 2. (II)

it is negative semi-definite; 3. (III)

it enjoys the following bounds:

[TABLE] 4. (IV)

it is homogeneous w.r.t. $Q$ with degree $-\frac{3}{2}$ and w.r.t. $S$ with degree $\frac{1}{2}$, i.e. $\,\bm{dT}_{aQ}^{S}=a^{-3/2}\bm{dT}_{Q}^{S}\,$ and $\,\bm{dT}_{Q}^{aS}=a^{1/2}\bm{dT}_{Q}^{S}\,$ for any $a>0$; 5. (V)

it is monotone w.r.t. $S^{1/2}QS^{1/2}$ (once range $S$ is fixed): $\,\bm{dT}_{Q_{0}}^{S_{0}}\preccurlyeq\bm{dT}_{Q_{1}}^{S_{1}}\,$ in the sense of self-adjoint operators on $\mathbb{H}(d)$ whenever $\,S_{0}^{1/2}Q_{0}S_{0}^{1/2}\preccurlyeq S_{1}^{1/2}Q_{1}S_{1}^{1/2}\,$ and $\,\operatorname{range}(S_{0})=\operatorname{range}(S_{1})\,$; in particular, $\bm{dT}_{Q}^{S}$ is monotone w.r.t. $Q\in\mathbb{H}_{++}(d)$ for fixed $S$.

Proof.

Slightly changing notations, we rewrite (A.1) as

[TABLE]

where matrices $\,U,U^{*}\,$ and $\,\varLambda\,$ come from Lemma A.2 and

[TABLE]

(I) Self-adjointness

Consider a scalar product

[TABLE]

We now introduce a following notation

[TABLE]

Then the above equality can be continued as follows:

[TABLE]

where $\,r=\operatorname{rank}(S)\,$. Thus the operator is self-adjoint.

(II) Boundedness and (III) eigenvalues

Denoting $\,\varDelta^{X}\,$ by $\,\varDelta\,$ (i.e. now $\,\varDelta=US^{1/2}XS^{1/2}U^{*}\,$) and taking into account the above expansion of an inner product, one obtains

[TABLE]

Note, that the function $\,f(\lambda_{i},\lambda_{j})\stackrel{{\scriptstyle\operatorname{def}}}{{=}}\frac{\sqrt{\lambda_{i}\lambda_{j}}}{\sqrt{\lambda_{i}}+\sqrt{\lambda_{j}}}\,$ is monotonously increasing in both arguments $\,\lambda_{i}\,$ and $\,\lambda_{j}\,$, thus

[TABLE]

For the sake of simplicity we introduce a new variable

[TABLE]

its Frobenius norm is written as

[TABLE]

Moreover, the following inequality for trace holds:

[TABLE]

Here $\,\varPi_{S}\,$ is the orthogonal projector onto range of $\,S\,$.

Then combining (A.2) with (A.3), the upper and lower bounds can be obtained as follows:

[TABLE]

Note, that if $\,S\,$ is degenerated, the lower bound becomes trivial.

(IV) Homogeneity and (V) monotonicity

Homogeneity follows directly from definition (A.1). Now we prove monotonicity. As range of $\,S^{1/2}QS^{1/2}\,$ is fixed, we may assume $\,S\succ 0\,$. Consider

[TABLE]

with replacement $\,M=S^{1/2}QS^{1/2}\,$ to be change of variables. As soon as $\,X\,$ is supposed to be fixed, it is enough to show that the differential $\,\bm{d_{M}}M^{-1/2}\,$ is monotone in $\,M\,$. Notice that the operator $\,\left(\bm{d_{M}}M^{-1/2}\right)^{-1}\,$ at point $M$ is equal to the differential of the inverse map $\,P\mapsto P^{-2}\,$ at point $\,P=M^{-1/2}\,$:

[TABLE]

In turn, $\,\bm{d_{P}}P^{-2}\,$ can be expressed as

[TABLE]

the right part of the above equation is self-adjoint, negative-definite and

[TABLE]

Choose $\,M_{1}\succcurlyeq M_{0}\succ 0\,$ (thus $\,M_{1}^{1/2}\succcurlyeq M_{0}^{1/2}\,$) and let $\,P_{i}=M_{i}^{-1/2}\,$ for $\,i=0,1\,$. Then for any fixed $\,X\in\mathbb{H}(d)\,$

[TABLE]

i.e. $\,\left.\bm{d_{P}}P^{-2}\right|_{P_{1}}\preccurlyeq\left.\bm{d_{P}}P^{-2}\right|_{P_{0}}\,$ and hence for the differential of $\,M\mapsto M^{-1/2}\,$ the inverse inequality holds: $\,\left.\bm{d_{M}}M^{-1/2}\right|_{M_{0}}\preccurlyeq\left.\bm{d_{M}}M^{-1/2}\right|_{M_{1}}\,$. This entails monotonicity of $\,\bm{dT}_{Q}^{S}\,$. ∎

Corollary A.1**.**

We define a following rescaled operator

[TABLE]

Then a following bound on its eigenvalues hold:

[TABLE]

Proof.

Notice that inequalities

[TABLE]

are a trivial consequence of Lemma A.3 (III). Now defining for any $\,1\leq k\leq\operatorname{rank}(S)\,$

[TABLE]

we obtain from (A.2) that

[TABLE]

Therefore, the above inequalities are sharp. ∎

Lemma A.4**.**

For any $\,Q_{0},Q_{1}\in\mathbb{H}_{++}(d)\,$, $\,S\in\mathbb{H}_{+}(d)\,$ consider

[TABLE]

Then

[TABLE]

Moreover, if $\,\norm{Q^{\prime}-I}<1\,$, then

[TABLE]

Remark 2**.**

The above inequality might seem confusing due to the fact that $\,\lambda_{min}(\cdot)\leq\lambda_{max}(\cdot)\,$, however this is explained by the fact that $\,\bm{dT}_{Q}^{S}\,$ is negative definite.

Proof.

Notice that

[TABLE]

Monotonicity and homogeneity with degree $\,-\frac{3}{2}\,$ of $\,\bm{dT}_{Q}^{S}(\cdot)\,$ (see Lemma A.3) yield

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

and respectively,

[TABLE]

The rest of the Lemma follow from the fact that

[TABLE]

and inequalities

[TABLE]

∎

A.2 Properties of $\,d_{BW}(Q,S)\,$

The next lemma ensures strict convexity of $\,d_{BW}(Q,S)\,$. In essence, the proof mainly relies on Theorem 7 Bhatia et al. (2018).

Lemma A.5**.**

For any $\,S\in\mathbb{H}_{+}(d)\,$ function $\,Q\mapsto d_{BW}^{2}(Q,S)\,$ is convex on $\,\mathbb{H}_{+}(d)\,$. Moreover, if $\,S\succ 0\,$, then it is strictly convex.

Proof.

According to (Bhatia et al., 2018, Theorem 7) function $\,h(X)=\operatorname{tr}X^{1/2}\,$ is strictly concave on $\,\mathbb{H}_{+}(d)\,$, hence function

[TABLE]

is convex on $\,\mathbb{H}_{+}(d)\,$ for any positive semi-definite $\,S\,$. Moreover, if $\,S\succ 0\,$, then $\,Q\mapsto S^{1/2}QS^{1/2}\,$ is an injective linear map, and therefore $\,d_{BW}^{2}(Q,S)\,$ is strictly convex. ∎

Further we present differentiability of $\,d_{BW}(Q,S)\,$ and its quadratic approximation.

Lemma A.6**.**

For any $\,Q\in\mathbb{H}_{++}(d)\,$, $\,S\in\mathbb{H}_{+}(d)\,$ function $\,d_{BW}^{2}(Q,S)\,$ is twice differentiable in $Q$ with

[TABLE]

Moreover, the following quadratic approximation holds: for any $\,Q_{0},Q_{1}\in\mathbb{H}_{++}(d)\,$

[TABLE]

with $\,Q^{\prime}\,$ defined in (A.5).

Proof.

Note that

[TABLE]

where $\,\delta\,$ comes from Lemma A.2. Furthermore, Lemma A.1 implies that

[TABLE]

Consequently, $\,d_{BW}^{2}(Q,S)\,$ is differentiable and

[TABLE]

Applying Lemma A.2 one obtains

[TABLE]

Quadratic approximation

Let $\,Q_{0},Q_{1}\in\mathbb{H}_{++}(d)\,$, $\,Q_{t}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}(1-t)Q_{0}+tQ_{1}\,$, $\,t\in[0,1]\,$. The Taylor expansion in integral form applied to $\,d_{BW}^{2}(Q_{t},S)\,$ implies

[TABLE]

Following the same ideas as in the proof of Lemma A.4 one obtains that

[TABLE]

and

[TABLE]

Thus

[TABLE]

∎

A.3 Central limit theorem for $\,Q_{n}\,$

First let us prove uniqueness and positive-definiteness of Bures–Wasserstein barycenter.

Proof of Theorem 2.1.

Since

[TABLE]

and $\,d_{BW}(Q,S)\to\infty\,$ as $\,\norm{Q}\to\infty\,$, the barycenter $\,Q_{*}\,$ always exists by continuity and compactness argument. In case $\,\mathbb{P}(\mathbb{H}_{++}(d))>0\,$ applying Lemma A.5 we obtain strict convexity of the integral

[TABLE]

and therefore, uniqueness of the minimizer $\,Q_{*}\,$.

To prove that $\,Q_{*}\succ 0\,$ consider arbitrary degenerated $\,Q_{0}\in\mathbb{H}_{+}(d)\cap\mathbb{A}\,$, $\,Q_{1}\in\mathbb{H}_{++}(d)\cap\mathbb{A}\,$ (which exists by Assumption 3) and $\,S\in\mathbb{H}_{++}(d)\,$. Let us define $\,Q_{t}=(1-t)Q_{0}+tQ_{1}\in\mathbb{A}\,$. We are going to show, that

[TABLE]

Consider eigen-decomposition $\,S^{1/2}Q_{0}S^{1/2}=U^{*}\varLambda U\,$, $\,\varLambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{r},0,\dots,0)\,$, where $\,r=\operatorname{rank}(Q_{0})\,$. Respectively, we write $\,C=US^{1/2}Q_{1}S^{1/2}U^{*}\,$ in a block form

[TABLE]

Cramer’s rule for inverse matrix and Laplace’s formula yield

[TABLE]

therefore

[TABLE]

Consequently,

[TABLE]

In the same way one can obtain

[TABLE]

Consequently,

[TABLE]

Since $\,\mathbb{P}(\mathbb{H}_{++}(d))>0\,$ by Assumption 3 and $\,d_{BW}^{2}(Q,S)\,$ is convex, we conclude that

[TABLE]

thus $\,Q_{0}\,$ cannot be a barycenter of $\mathbb{P}$. This yields $\,Q_{*}\succ 0\,$.

Since $\mathcal{V}(\cdot)$ is convex and barycenter of $\mathbb{P}$ is positive-definite and unique, it is characterized as a stationary point of Fréchet variation on subspace $\mathbb{A}$, i.e. as a solution to equation

[TABLE]

as required. ∎

The proof of CLT widely uses the concept of covariance operators on the space of optimal transportation maps and on the space of covariance matrices.

Covariance operator on the space of optimal maps

Consider $\,T_{i}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}T_{Q_{*}}^{S_{i}}\,$ and $\,T^{n}_{i}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}T_{Q_{n}}^{S_{i}}\,$. We define covariance $\,\bm{\varSigma}\,$ of $\,T_{i}\,$, its empirical counterpart $\,\bm{\varSigma_{n}}\,$, and its data-driven estimator $\,\bm{\hat{\varSigma}_{n}}\,$ as follows

[TABLE]

Covariance operators on the space of covariance matrices

Let $\,Q_{n}\,$ be an empirical barycenter. The covariance of $\,Q_{n}\,$ and its empirical counterpart are defined as

[TABLE]

[TABLE]

where

[TABLE]

Now we are ready to prove the central limit theorem for the empirical barycenter $\,Q_{n}\,$ (Theorem 2.2). However, for the sake of transparency we provide below a complete statement.

Theorem** (Central limit theorem for the covariance of empirical barycenter).**

The approximation error rate of the Fréchet mean $\,Q_{*}\,$ by its empirical counterpart $\,Q_{n}\,$ is

[TABLE]

Moreover, if $\,(\bm{\varSigma})_{\mathbb{M}}\,$ is non-degenerated, then

[TABLE]

Proof of Theorem 2.2.

The proof consists of two parts: proof of (A) and (B).

Proof of (A)

As $\,\mathcal{V}_{n}(\cdot)\,$ are convex functions, they a.s. uniformly converge to strictly convex function $\,\mathcal{V}(\cdot)\,$ on any compact set by the uniform law of large numbers. Therefore, their minimizers also converge a.s. $\,Q_{n}\to Q_{*}\,$. In particular, $\,Q_{n}\succ 0\,$ with dominating probability.

Applying the expansion from Lemma A.2 at point $\,Q_{*}\,$ implies

[TABLE]

where $\,Q_{t}=(1-t)Q_{*}+tQ_{n}\,$ and $\,\alpha(\cdot;Q_{n},S_{i})\,$ is a self-adjoint operator on $\,\mathbb{H}(d)\,$ s.t. \,\left\langle X,\alpha\bigl{(}X;Q_{n},S\bigr{)}\right\rangle=o\bigl{(}\bigl{\langle}X,\bm{dT}_{Q_{*}}^{S}(X)\bigr{\rangle}\bigr{)}\, as $\,Q_{n}^{\prime}\to I\,$ uniformly in $\,S\,$ and $\,X\,$ according to Lemma A.4. Note, that the condition for $\,Q_{n}\,$ being a barycenter is $\,\bm{\varPi}_{\mathbb{M}}\left(\frac{1}{n}\sum_{i}T^{n}_{i}-I\right)=0\,$. This fact together with averaging of (A.11) over $\,i\,$ give:

[TABLE]

where

[TABLE]

and $\,\bm{F}_{n}\,$ is defined in (A.9). Recall that $\,\bm{F}\,$ (A.9) is a population counterpart of $\,\bm{F_{n}}\,$. This operator is correctly defined since by Lemma A.3 one can show that it is self-adjoint, positive definite and bounded:

[TABLE]

Since by the law of large numbers $\,\bm{F}_{n}\to\bm{F}\,$ and \,\alpha_{n}(X)=o\bigl{(}\norm{\bm{F}_{n}(X)}\bigr{)}\,, we obtain from (A.12)

[TABLE]

where $\,\bm{F}^{-1}\,$ is a bounded linear operator, because $\,\bm{dT}_{Q_{*}}^{S}\,$ is negative definite for any $\,S\succ 0\,$ by Lemma A.3. The result (A) follows immediately from the CLT for $\,\bm{\varPi}_{\mathbb{M}}\overline{T}_{n}\,$.

Proof of (B)

Note that result (A) is equivalent to the fact, that

[TABLE]

To ensure convergence of $\,\bm{\hat{\Upxi}_{n}}\rightarrow\bm{\Upxi}\,$ we need to show that

• •

$\,\bm{\hat{\varSigma}_{n}}\rightarrow\bm{\varSigma}\,$ (is proved in Lemma B.2);

• •

$\,\bm{\hat{F}_{n}}\rightarrow\bm{F}\,$.

Convergence of $\,\bm{\hat{F}_{n}}\,$ to $\,\bm{F}\,$

Monotonicity and homogeneity with degree $\,-\frac{3}{2}\,$ of $\,\bm{dT}_{Q}^{S}\,$ (see Lemma A.3) yield

[TABLE]

where $\,Q^{\prime}_{n}\,$ comes from (A.5). This naturally leads to the following relation:

[TABLE]

Note, that $\,\lambda_{min}(Q^{\prime}_{n})\geq 1-\norm{Q^{\prime}_{n}-I}\,$ and $\,\lambda_{max}(Q^{\prime}_{n})\leq 1+\norm{Q^{\prime}_{n}-I}\,$. This immediately implies

[TABLE]

Since $\,\bm{F_{n}}\to\bm{F}\,$ this implies $\,\bm{\hat{F}_{n}}\to\bm{F}\,$.

The above results ensures validity of substitution $\,\bm{\Upxi}\,$ by $\,\bm{\hat{\Upxi}_{n}}\,$. This yields (B). ∎

The asymptotic convergence results for $\,d_{BW}(Q_{n},Q_{*})\,$ is a straightforward corollary of the above theorem. Here is the proof.

Proof of Corollary 2.1.

Since $\,Q_{n}\to Q_{*}\,$ a.s., Lemma A.6 implies

[TABLE]

By Theorem 2.2 $\,\sqrt{n}(Q_{n}-Q_{*})\,$ is asymptotically normal and centred, therefore

[TABLE]

where $\,Z\in\mathbb{M}\subset\mathbb{H}(d)\,$ and $\,Z\sim\mathscr{N}\left(0,\bm{\Upxi}\right)\,$.

Note, that $\,Q_{n}\rightarrow Q_{*}\,$, $\,\bm{\hat{\Upxi}}_{n}\rightarrow\bm{\Upxi}\,$, and $\,\bm{dT}_{Q_{n}}^{Q_{n}}\rightarrow\bm{dT}_{Q_{*}}^{Q_{*}}\,$. The latter result follows from Lemma A.3 (IV, V):

[TABLE]

This yields

[TABLE]

where $\,Z_{n}\sim\mathscr{N}\left(0,\bm{\hat{\Upxi}}_{n}\right)\,$. This, in turn, entails

[TABLE]

Appendix B Concentration of barycenters

The next lemma is a key ingredient in the proof of concentration result for $\,Q_{n}\,$.

Lemma B.1**.**

Consider

[TABLE]

where

[TABLE]

Then

[TABLE]

whenever $\,\eta_{n}<\frac{4}{3}\,$ and $\,Q_{n}\succ 0\,$.

Proof.

Let us define $\,Q_{t}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}tQ_{n}+(1-t)Q_{*}\,$ for $\,t\in[0,1]\,$. Due to Lemmas A.3 and A.4 we have for any $\,S\in\mathbb{H}_{+}(d)\,$

[TABLE]

Therefore,

[TABLE]

At the same time,

[TABLE]

Hence

[TABLE]

Rewriting the inequality above we obtain

[TABLE]

provided that $\,\eta_{n}<\frac{4}{3}\,$. ∎

Proposition 1** (Concentration of $\,\overline{T}_{n}\,$; Hsu et al. (2012), Theorem 1).**

Under Assumption 2 it holds

[TABLE]

Before proving concentration results, we define operator $\,\bm{F^{\prime}}(X)\,$ as:

[TABLE]

Proof of Theorem 2.3.

Let $\,t_{n}\,$ be s.t. the following upper bound on $\,\gamma_{n}(t_{n})\,$ from Lemma B.3 holds:

[TABLE]

It is easy to see that this condition is fulfilled for $t_{n}=nt_{F}-\log(m)$ under a proper choice of generic constant in definition of $t_{F}$. Then with $\,\mathbb{P}\geq 1-2me^{-nt_{F}}\,$ a following bound holds

[TABLE]

The above facts together with definition of $\,\eta_{n}\,$ B.1 yield

[TABLE]

Combining the above bounds with Proposition 1, we obtain:

[TABLE]

Now it follows from Lemma B.1 that

[TABLE]

whenever $\,\frac{c_{Q}}{2\sqrt{n}}(d+t)\leq\frac{2}{3}\,$. Here we used that $\,Q_{n}\succ 0\,$ if at least one of matrices $\,S_{1},\dots,S_{n}\,$ is non-degenerated. ∎

Proof of Corollary 2.2.

To prove this result we use Lemma A.6 and choose $\,Q_{0}=S=Q_{*}\,$, $\,Q_{1}=Q_{n}\,$. Thus we obtain

[TABLE]

Hence by Theorem 2.3

[TABLE]

with probability at least $\,1-2me^{-nt_{F}}-e^{-t^{2}/2}-(1-p)^{n}\,$. ∎

B.1 Central limit theorem and concentration for $\,\mathcal{V}_{n}\,$

Proof of Theorem 2.4.

By definition empirical Fréchet variance is

[TABLE]

Lemma A.6 ensures the following bound on $\,\mathcal{V}_{n}(Q_{*})-\mathcal{V}_{n}(Q_{n})\,$:

[TABLE]

with $\,Q^{\prime}_{n}\,$ defined in (2.4). The above quadratic bound together with $\,Q_{n}\to Q_{*}\,$, $\,\bm{F}_{n}\to\bm{F}\,$ and $\,\sqrt{n}(Q_{n}-Q_{*})\rightharpoonup\mathscr{N}(0,\bm{\Upxi})\,$ yield:

[TABLE]

On the other hand, by classical central limit theorem we obtain:

[TABLE]

∎

Proof of Theorem 2.5.

Following the proof of Theorem 2.4 we consider $\,\mathcal{V}_{n}(Q_{*})-\mathcal{V}_{n}(Q_{n})\,$:

[TABLE]

Following the proof of Theorem 2.3, we obtain that with $\,\mathbb{P}\geq 1-2me^{-t_{F}n}-e^{-t^{2}/2}-(1-p)^{n}\,$ the following upper bounds hold:

[TABLE]

Thus

[TABLE]

and consequently

[TABLE]

Now we consider a difference $\,\mathcal{V}_{n}(Q_{*})-\mathcal{V}(Q_{*})\,$. According to Assumption 2 $\,S\,$ and therefore $\,d_{BW}^{2}(Q_{*},S)\,$ are sub-exponential r.v. with some parameters $\,(\nu,b)\,$. Then Lemma B.4 ensures

[TABLE]

with probability $\,1-2e^{-t^{\prime}}\,$. Combining two above bounds, we obtain:

[TABLE]

with probability

[TABLE]

Choosing $\,t^{\prime}=t^{2}/2\,$, we obtain

[TABLE]

∎

B.2 Auxiliary results

Lemma B.2**.**

Let $\,\norm{Q_{n}^{\prime}-I}\leq\frac{1}{2}\,$; then

[TABLE]

where

[TABLE]

where $\,\kappa(Q_{*})={\norm{Q_{*}}}{\norm{Q_{*}^{-1}}}\,$ is the condition number of matrix $\,Q_{*}\,$ and $\,\norm{\cdot}_{1}\,$ is $1$-Schatten (nuclear) norm of an operator.

Proof.

Note, that for any $\,\left(T^{n}_{i}-I\right)\otimes\left(T^{n}_{i}-I\right)\,$ the following decomposition holds

[TABLE]

Summing by $\,i\,$ yields

[TABLE]

Note, that each

[TABLE]

Lemmas A.3 (III) and A.4 yield

[TABLE]

hence $\,\frac{1}{n}\sum_{i}\norm{T^{n}_{i}-T_{i}}_{F}^{2}\leq\beta_{n}^{2}\,$. The above expression together with (B.6) and Cauchy-Schwarz inequality lead to the upper bound on $\,\norm{\bm{\hat{\varSigma}_{n}}-\bm{\varSigma_{n}}}_{1}\,$:

[TABLE]

∎

Further we present concentration of $\bm{F}_{n}$ around $\bm{F}$. Denote as $\,\norm{\cdot}_{\psi_{2}}\,$ an Orlicz norm with Young function $\,\psi_{2}(x)=e^{x^{2}}-1\,$, i.e.

[TABLE]

Then sub-Gaussianity of a r.v. $\,X\,$ is equivalent to $\,\norm{X}_{\psi_{2}}<\infty\,$ and it ensures

[TABLE]

Lemma B.3** (Concentration of $\bm{F^{\prime}}_{n}$, Proposition 2 in Koltchinskii et al. (2011)).**

There exists a constant $\,\mathtt{C}>0\,$, s.t. for all $\,t>0\,$ it holds with probability at least $\,1-e^{-t}\,$

[TABLE]

where $\,\sigma_{F}\stackrel{{\scriptstyle\operatorname{def}}}{{=}}\norm{\mathbb{E}\left(\bm{dt}_{Q_{*}}^{S}-\bm{F^{\prime}}\right)^{2}}\,$, $\,U\stackrel{{\scriptstyle\operatorname{def}}}{{=}}\norm{\norm{\bm{dt}_{Q_{*}}^{S}-\bm{F^{\prime}}}}_{\psi_{2}}\,$.

Lemma B.4** (Sub-exponential tail bounds).**

Suppose that $\,X\,$ is sub-exponential with parameters $\,\nu,b\,$. Then

[TABLE]