Gaussian approximation for
empirical barycenters

Nazar Buzun

<EMAIL_ADDRESS>Institute of Science and Technology, Bolshoy Boulevard 30, bld. 1
Moscow, Russia 121205

In this work we consider  Wasserstein barycenters (average in Wasserstein distance) in Fourier basis. We prove that random Fourier parameters of the barycenter converge to Gaussian random vector by distribution. The convergence approximation has been done in finite-sample condition with  rate $O(p/\sqrt{n})$ depending on measures count ($n$) and the dimension of parameters ($p$).

Keywords: GAR, Wasserstein distance, multivariate central limit theorem, statistical learning, convex analysis.

§ INTRODUCTION

Monge-Kantorovich distance or Wasserstein distance is a distance between measures. It represents a transportation cost of measure $\mu_1$ into the other measure $\mu_2$.

\[
W_p(\mu_1, \mu_2) =  \min_{\pi \in \Pi[\mu_1, \mu_2] } \left(  \int \| x - y \|^p d \pi (x, y) \right)^{1/p}
\]

where the condition $\pi \in \Pi[\mu_1, \mu_2]$ means that $\pi(x,y)$ has two marginal distributions:  $\int_y d \pi(x,y) = d\mu_1(x)$ and $\int_x d \pi(x,y) = d \mu_2(y)$. We focus on regularized W1 distance with probabilistic space $\{\R^d, \mathcal{B}(\| \cdot \|_2), L^1\}$

\[
\widetilde{W}_1(\mu_1, \mu_2) =  \min_{\pi \in \Pi[\mu_1, \mu_2] }  \int \| x - y \| d \pi (x, y)  + R_{\varepsilon}(\pi)
\]
where $R_{\varepsilon}(\pi)$ is a relatively small addition which improves differential properties of the distance. Namely without $R_{\varepsilon}(\pi)$ we can only bound the first derivative, with it we can bound also the second derivative.
There is the notion of mean in Wasserstein distance, called barycenter $\muhat$. And it is the main object in this paper.

\[
\muhat = \argmin_{\mu} \sum_{i=1}^n \widetilde{W}_1(\mu, \mu_i)
\]

Barycenters are center-of-mass generalization. If we look at the barycenter of a set of uniform measures it fits the common structure form of these measures. If the measures are sampled from some distribution then their barycenter can be treated as an empirical approximation of the distribution mean. A simple example is a circles set with means $\{m_i\}$ and radius's $\{r_i\}$.

Illustration for $W_2$ distance computation between two circles $(m_1, r_1)$ and $(m_2, r_2)$.

\[
W_2^2 \big((m_1, r_1), (m_2, r_2) \big) = \frac{1}{2\pi} \int_{0}^{2\pi}  [(m_2 - m_1) - (r_2 - r_1) \cos(a)]^2 + [(r_2 - r_1) \sin(a)]^2 da
\]

\[
= (m_2 - m_1)^2 + (r_2 - r_1) ^2
\]

Their $W_2^2$ barycenter is also c circle with mean $m = \frac{1}{n} \sum_{i=1}^n m_i$ and radius $r = \frac{1}{n} \sum_{i=1}^n r_i$. We  refer to papers [1], [7] for the overview of the barycenters and related study.

It is well known that the center-of-mass in $l_2$ norm converges to a Gaussian random vector. As for the barycenter it is also expected to have some Gaussian properties. For example the measures are Gaussian themselves or one-dimensional or circles set then the Gaussian structure of the barycenter is evident. In circles set case the mean and radius converges to some Gaussian variables as a sum of independent observations according to Central Limit Theorem.  In one-dimensional case denoting distribution functions by $F_i(x)$

\[
W_1(\mu_1, \mu_2)  =  \int_0^1 | F_1^{-1}(s) - F_2^{-1}(s) | ds
\]

one get

\[
\widehat{F}^{-1}(s) = \frac{1}{n} \sum_{i=1}^n F_i^{-1}(s)
\]

In the case with Gaussian measures with zero mean and variances $\{S_i\}$

\[
W_2^2(\mu_1, \mu_2) =  \tr(S_1) + \tr(S_2) - 2 \tr((S_2^{1/2}S_1S_2^{1/2})^{1/2})
\]

and for some non-random matrix $S_*$ [11]

\[
\widehat{S} = \frac{1}{n} \sum_{i=1}^n (S_*^{1/2}S_iS_*^{1/2})^{1/2} + O(1/n)
\]

In both last cases one have a mean of independent random variables that converges to Gaussian random variable (or to Gaussian process in case of $\widehat{F}^{-1}(s)$ by Donsker's Theorem).
In general case it appears to be very difficult  to reveal such convergence because the barycenter doesn't have an explicit equation and it is an infinite-dimensional object. In order  to handle with this difficulty  we propose to construct a sum of independent variables using projection into Fourier basis and some novel results from statistical learning theory. The perspective of Fourier Analysis provides a suitable representation of the Wasserstein distance and it is already studied in the literature [10]. Denote a range of size $p$ of the barycenter Fourier coefficients by

\[
\opttheta =  \mathcal{F}_p \left( \frac{d\muhat(x)}{dx} \right)
\]

The first our result is that for some non-random matrix $D$, non-random vector $\thetas$ and independent random vectors $\{ \xi_i \}$

\[
\left\| D \left( \opttheta  - \thetas \right) - \sum_{i=1}^n \xi_i \right\|  = O\left( \frac{p}{\sqrt{n}} \right)
\]

Further we show that for some Gaussian vector $Z$

\[
W_1 \left( D ( \opttheta  - \thetas ), Z \right ) = O\left( \frac{p^{3/2}}{\sqrt{n}} \right)
\]

\[
\left | P \left( \| D ( \opttheta  - \thetas ) \| > x \right) -  P \left(\|Z\| > x \right) \right | = O\left( \frac{p}{\sqrt{n}} \right)
\]

Statistical Application: The last statement allows as to obtain the confidence region of parameter $\opttheta$ and describe the distribution inside the region. The bootstrap procedure  validity
follows from our proof as well. If one sample $\| D ( \theta^{boot} -  \opttheta ) \|$ using bootstrap it would be close by quantiles to the random variable $\| D ( \opttheta - \thetas ) \|$. This is also relates to the construction of the confidence region.

The Structure of this paper is following. The main Theorems are in Section 2. In Sections 5,6 we compute derivatives of the Wasserstein distance using infimal convolution of support functions. Section 3 deals with independent parametric models and describes how one can approximate parameter deviations by a sum of independent random vectors $\{ \xi_i \}$. In Section 4 we explore the barycenters model and check the required assumptions from the 3-rd Section. The final part, Gaussian approximation of the parameter $\opttheta_p$, is completed in Section 7, where we prove that $\{ \xi_i \}$ is close to $Z$ by distribution and by $W_1$.

§ THE MAIN RESULT

Consider a set of random measures (random measure is a measure-valued random element) with densities $\phi_1, \ldots, \phi_n$.  Let the barycenter measure $\muhat$ has density $\widehat{\phi}$ and Fourier coefficients $\opttheta = \theta(\widehat{\phi}) \in \R^{\infty}$.

\[
\widehat{\phi}=\argmin_{\phi}\sum_{i=1}^{n} \widetilde{W}_1(\phi, \phi_i)
\]

Let Fourier basis $\{ \psi_k \}_{k=1}^{\infty}$ has a Gram function of the scalar product $G(x)$, such that for some function $f$

\[
\langle f,\psi_{k}\rangle_{G}=\int f(x)\psi_{k}(x) G(x)dx
\]

\[
\theta(\varphi)[k]=\int\varphi(x)\psi_{k}(x)dx
\]

Denote Fourier coefficients of the other measures $\forall i: \theta_i = \theta(\varphi_i) \in \R^{\infty}$. Basing on Lemma <ref> define an independent parametric model with dataset ($\theta_1, \ldots, \theta_n$) and parameter $\theta$.

\[
L(\theta) = \sum_{i=1}^n l(\theta - \theta_i)
\]

\[
l(\theta - \theta_i) = \max_{\eta \in \bigcap \mathcal{E}_{x}}\langle\eta, \theta - \theta_i\rangle - \varepsilon \eta^T (K \circ G) \eta =
\widetilde{W}_1(\phi, \phi_i)
\]

and  $\bigcap \mathcal{E}_{x}$ is a Sobolev ellipsoids intersection. Each ellipsoid $\mathcal{E}_x$ has matrix $K_x = \nabla^T \psi \nabla \psi^T (x)$ such that

\[
\left( \sum_{k \in \N^{d}_{+} }\eta_{k} \nabla \psi_k(x) \right)^{2}=\eta^{T} K_x \eta
\]

\[
\bigcap \mathcal{E}_{x}=\left\{ \eta : \forall x: \eta^{T} K_x \eta \leq 1 \right\}
\]

Define a positive matrix $K \circ G = \int K_x G(x) dx$ such that in case $\psi_{k}(x) = e^{i k^T x / T}$

\[
K \circ G =   \left(\begin{matrix}
1/T^2 & 0 & 0 &  \ldots 0 \\
0  &  \ddots & 0 &  \ldots 0 \\
0 & \ldots &  k^2/T^2  & \ldots 0 \\
\ldots & \ldots & \ldots & \ldots
\end{matrix} \right)
\]

Define for this model MLE parameter value and reference parameter value:

\[
\widehat{\theta}=\argmin_{\theta}L(\theta)
\]

\[
\theta^{*}=\argmin_{\theta}\E L(\theta)
\]

Define a local region around $\thetas$

\[
\Omega(\rr) = \{ \theta: \|D (\theta - \thetas) \| \leq \rr \}
\]

where $D$ is the Fisher matrix of the model

\[
D^2 = - \nabla^2 \E L(\thetas)
\]

Let the random Fourier parameters of the dataset have a common density $\theta_1 \ldots \theta_n \sim q(\theta)$ and it fulfills condition $\forall \theta \in \Omega(\rr)$

\[
\int   \| D^{-1} \nabla q(\theta) \| d\theta  = \frac{C_Q}{\sqrt{n}}
\]

Let $\opttheta, \thetas \in \R^{\infty}$ are Fourier coefficients of the MLE and reference barycenter defined above, then with probability $\hiprobt$

\[
\normp{ D (\opttheta - \theta^*) -  D^{-1} \nabla L (\thetas) }
\leq  \rombt
\]

where $\rombt$ is defined in Section <ref> and has asymptotic

\[
\rombt  =  \frac{ \sqrt{n} O (\rr C_Q + \rr \sqrt{p_D} + \sqrt{2t}) }{\varepsilon \lambda_{\min} \big(D K \circ G D \big)}
+ o \left( \frac{   1  }{\sqrt{n}} \right)
\]

and $p_D$ is an ellipsoid entropy (Section <ref>) with matrix $D$

\[
p_D = \sqrt{\sum_i \frac{ \log^2 (\lambda_i^2(D)) }{\lambda_i^{2}(D)}  }
\]

and with probability $\hiprobt$

\[
\rr \leq 4 \|D^{-1} \nabla L(\thetas) \| \leq  \frac{ 8 \sqrt{n} (1 + \sqrt{2t}) }{\lambda_{\min}^{1/2} \big( D K\circ G D \big)}
\]

Basing on Theorem <ref> one have to prove Assumptions 1,2,3 from which follows

\[
\normp{ D (\opttheta - \theta^*) -  D^{-1} \nabla L (\thetas) }
\leq \{\delta(\rr)  +  \zz(t)\} \rr = \rombt
\]

with probability $\hiprobt$. The Assumptions 1,2,3 are proven in Section <ref> where also is shown that

\[
\delta(\rr) = \frac{ \rr n }{\varepsilon \lambda_{\min} \big(D K \circ G D \big)} \int  \| D^{-1} \nabla q(\theta) \| d\theta
\]

\[
\zz(t) =  E +  \sqrt{2t(\textbf{v}^{2}+2R E)}   +\frac{tR}{3}
\]

\[
E = 12 \mathbf{v} \sqrt{ 2 p_D} + 24 R    p_D
\]

\[
\textbf{v}^2  = \frac{  n }{ \varepsilon^2\lambda_{\min}^2 \big(D K \circ G D \big)}
\]

\[
R = \frac{ 1 }{\varepsilon \lambda_{\min} \big(D K \circ G D \big)}
\]

Setting $\textbf{v}$ and $R$ in previous equations gives an asymptotic

\[
\zz(t) = \frac{\sqrt{n} (12 \sqrt{2p_D} + \sqrt{2t} ) }{ \varepsilon \lambda_{\min} \big(D K \circ G D \big) }
+ O \left( \frac{1}{n} \right)
\]
Lemma <ref> gives bound

\[
\| D^{-1} \nabla l \| \leq \frac{1}{\lambda_{\min}^{1/2} \big( D  K\circ G D \big)}
\]

From this bound Hoefding's inequality [5] follows bound for $\| D^{-1} \nabla L(\thetas) \|$.

Define additional Fisher matrix corresponded to the projection into first $p$ elements of the parameter $\theta$ (ref. for details in Section <ref>).

\[
{\breve{D}}^2  =  D^2_{p\times p} - D^2_{p\times \infty} D^{-2}_{\infty\times \infty} D^2_{\infty\times p}
\]

such that

\[
D^2 =  \left(
\begin{matrix}
D^2_{p \times p} & D^2_{p \times \infty} \\
D^2_{\infty \times p} & D^2_{\infty \times \infty}
\end{matrix} \right)
\]

and define the gradient of the projection into first $p$ elements of the parameter $\theta$.

\[
\scorer = \nabla_{1 \ldots p} -   D^2_{p\times \infty} D^{-2}_{\infty\times \infty} \nabla_{p \ldots \infty}
\]

Let $\opttheta_p, \thetas_p \in \R^p$ are  the first $p$ Fourier coefficients of the MLE and reference barycenters,  and $Z$ is a Gaussian vector $\ND(0,  \Var [{\breve{D}}^{-1} \scorer L(\thetas)])$. Then with probability $(\hiprobt)$ $W_1$ and probability distances to $Z$ are bounded as follows

\[
W_1 ({\breve{D}} (\opttheta_p - \thetas_p), Z) \leq  \mu_3 \,    O(\log(n)) + \rombt
\]
and  $\forall z \in \R_{+}$

\[
| \P ( \| D (\opttheta - \thetas) \| > z )  -  \P ( \| Z \| > z ) |
\leq  C_A \big(\mu_3 \,    O(\log^2 n) + \rombt \big)
\]

where $\rombt$ is defined in Theorem <ref>,  $C_A = O(1/z)$  is anti-concentration constant defined in Section <ref> and

\[
\mu_3 \leq  \frac{ 4\sqrt{2} p   }{\lambda^{1/2}_{\min} \big( D  K\circ G D \big)}
\]

Bind Theorems <ref> and <ref>.  Form Theorem <ref> follows that the bound in Theorem <ref> also holds for projection of the parameter $\theta$:

\[
\|  \breve{D} (\opttheta_p - \thetas_p) -  \breve{D}^{-1}\scorer L(\thetas) \| \leq \rombt
\]

So  with probability $\hiprobt$

\[
W_1 ( \breve{D} (\opttheta_p - \thetas_p),  Z) = \min_{\pi(\opttheta, Z)} \E \| \breve{D} (\opttheta_p - \thetas_p) -  Z \|
\]

\[
\leq W_1 (\breve{D}^{-1}\scorer L(\thetas),  Z)   + \rombt
\]

Furthermore  from  Theorem <ref>  follows

\[
W_1(\breve{D}^{-1}\scorer L(\thetas), Z) \leq \sqrt{2} \mu_3 \left (1 + \log(2 \sqrt{ \tr \{ \Sigma \} } \mu_2) - \log (\mu_3) \right)
\]

where $\Sigma = \Var [\breve{D}^{-1}\scorer L(\thetas)]$ and setting $X_i =\breve{D}^{-1}\scorer l(\thetas - \theta_i) $

\[
\mu_3 =  \sum_{i=1}^n  \E \| \Sigma^{-1/2} (X_i - X'_i) \| \| \Sigma^{-1/2} X_i \| \| X_i - X'_i \| \leq
4  \max \| X_i \| \sum_{i=1}^n  \E X_i^T \Sigma^{-1} X_i
\]

\[
\sum_{i=1}^n  \E X_i^T \Sigma X_i     =  \tr \left\{ \Sigma^{-1} \sum_{i=1}^n  \E   X_i X_i^T \right\} = p
\]

\[
\max \| X_i \| = \| \breve{D}^{-1}\scorer l(\thetas - \theta_i)  \| \leq \| D^{-1}\nabla l(\thetas - \theta_i)  \|
\leq \frac{1}{ \lambda^{1/2}_{\min} \big( D  K\circ G D \big) }
\]

Analogically one can make a consequence from Theorems <ref> and <ref>.
Let $C_A$ is the anti-concentration constant of the distribution $\P(\| Z \| > x)$, then

| (  D̆ (_p - _p)  >  z )  -  (  Z  > z ) |

≤| (  D̆^-1 L()  >  z )  -  (  Z  > z  ) |
+ C_A

\[
| \P ( \| \breve{D}^{-1} \scorer L(\thetas) \| > z  )  -  \P ( \| Z \| >  z  ) | \leq  C_A \mu_3 O( \log^2n )
\]

As for the anti-concentration constant it can be estimated from Pinsker's inequality applied in the next equation

\[
\P(\| Z \|  \in [z, z + \Delta])  = \frac{1}{(2\pi \det \Sigma)^{1/2}} \int_{z < \| x \| < z + \Delta} e^{-x^T \Sigma^{-1} x/2} dx  = O\left( \frac{\Delta}{z}  \log(\tr\{\Sigma\}) \right)
\]

§ STATISTICAL LEARNING THEORY

In this section we consider an infinite dimensional statistical model $L (\theta) $. Let parameter $\theta $ consists of two parts $(u, v) $, such that $u=\theta_{1\ldots p} \in \R^p$. Given a finite dataset we are going to find MLE deviations basing on three assumptions listed below. Further we will specify these assumptions for independent models and apply to barycenter model from the previous section.

§.§ General approach

Let the Likelihood function $L(\theta) = L(\theta,Y)$  depends on parameters vector $\theta = (u,v)  $ and  a fixed dataset $Y$ of size $n$. Denote the parameter MLE and refernce values:

\[
\opttheta= \argmax_\theta L(\theta )
\]

\[
\thetas = \argmax_\theta \E L(\theta)
\]

We are going to study deviations of $\opttheta$ and $u$ in the following sense. For some  matrix $D$ and random vector $\xi$

* $ \| \opttheta - \thetas \| $ is expected to be of order $O(1/\sqrt{n})$.

* $  D(\opttheta - \thetas)  \approx \xi $

* $ L(\opttheta) - L(\thetas)  \approx \| \xi \|^2 / 2 $

Denote the stochastic part of the Likelihood

\[
\zeta(\theta) = L(\theta) - \E L(\theta)
\]

Involve the Fisher matrix

\[
D^2 = - \nabla^2 \E L(\thetas) =  \left(
\begin{matrix}
\Ftt & \Fte \\
\Fet & \Fee
\end{matrix} \right)
\]

It would be easier to deal with the model  if matrix $\Fbb$ has block-diagonal view ($\Fte = 0$). One can make parameter replacement in order to satisfy to this condition.  Define a new variable $\vartheta = \vartheta(u,  v)$ such that

\[
\nabla_\vartheta  \nabla_u^T \E L(\thetas) = \nabla_\vartheta  \nabla_u^T \E L(\thetas) = 0
\]

\[
\vartheta = v +  D^{-2}_{v} \Fet u,
\]

or in other words the parameters transformation matrix is

\[
S =
\left(\begin{matrix}
I & 0 \\
D^{-2}_{v} \Fet & I
\end{matrix} \right),
\quad
S^{-1} =
\left(\begin{matrix}
I &  0  \\
- D^{-2}_{v}\Fet  & I
\end{matrix} \right)
\]

The gradient in the new coordinates $(u, \vartheta)$  may be obtained by rule $\nabla(u, \vartheta)= (S^{-1})^T \nabla(u,  v)$. Use notation $\scorer$ for the first part of it

\[
\scorer = \nabla_{u} (u, \vartheta) = \nabla_{u} -   \Fte  D_{v}^{-2} \nabla_{v}
\]

The Fisher matrix after parameters replacement changes by rule $\Fbb(u, \vartheta) =  (S^{-1})^T \Fbb S^{-1}$,  so in the new coordinates it has view

\[
D^2(u, \vartheta) = -\nabla^2 \E L(u^*, \vartheta^*) =  \left(
\begin{matrix}
\breve{D}^2 & 0 \\
0 & D_{\vartheta}^2
\end{matrix} \right)
\]

\[
\breve{D}^2  = \Ftt - \Fte D_v^{-2} \Fet
\]
Define a local region around point $\thetas$

\[
\Omega(\rr) = \{\theta: \| D (\theta - \thetas) \| \leq \rr \}
\]

Now we write down three conditions on the Likelihood derivatives essential for the deviations of $\opttheta$. The first and second conditions holds in the local region $\Omega(\rr)$. The third condition is required to make expansion of local statements to the whole parameter space $\R^{\infty}$. Further we will show that from these conditions also follows deviation bounds of the parameter $\widehat{u}$ or  in other words from  deviations bound of $\opttheta$ follows bound of $\widehat{u}$.

Assumption 1: In the region $\Omega(\rr)$

- D^-1  { ∇L(θ) -  ∇L() } - D (θ- )      ≤    δ()
Assumption 2: In the region $\Omega(\rr)$ with probability $\hiprobt$

sup_θ∈Ω()   D^-1  { ∇ζ(θ) - ∇ζ() }       ≤     (t)
Assumption 3:
The Likelihood function is convex  $(-\nabla^2 L(v) \geq 0)$ or the expectation of Likelihood function is upper-bounded by a strongly convex function $(\E L(\vs) - \E L(v) \geq b \rr_v^2$).

Use notation

\[
\{\delta(\rr)  +  \zz(t)\} \rr = \rombt
\]

Let the Likelihood function is convex $(-\nabla^2 L(v) \geq 0)$ and for $\rr_v$ (assigned further) $\delta(\rr)  +  \zz(t) \leq 1/2$. Then under Assumptions 1, 2 with probability $\hiprobt$

\[
\rr \leq 4 \|D^{-1} \nabla L(\thetas) \|
\]

\[
\| D (\opttheta - \theta^*) -  D^{-1} \nabla L (\thetas) \|
\leq  \rombt
\]

\[
\| \breve{D} (\widehat{u} - u^*) -  \breve{D}^{-1} \scorer L (\thetas) \|
\leq \rombt
\]

From $(-\nabla^2 L(\theta) \geq 0)$ and ($L(\opttheta) > L(\thetas)$) follows that the local region $\Omega(\rr)$ that includes $\opttheta$ should covers the next region

\[
\Omega(\rr)  \supset \{\theta: L(\theta) \geq L(\thetas) \}
\]

Estimate the minimum possible radius of $\Omega(\rr)$ that satisfy the previous condition.

\[
0 \geq L(\thetas) - L(\theta) = -(\theta-\thetas)^T \nabla L(\thetas) - \frac{1}{2} (\theta - \thetas)^T \nabla^2 \zeta(\theta_0) (\theta - \thetas)  + \frac{1}{2} \| D(\theta_0) (\theta - \thetas) \|^2
\]

\[
\{ \text{with probability } (\hiprobt) \} \geq - \|D^{-1} \nabla L(\thetas) \| \rr - \frac{ \zz(t) }{2} \rr^2 + \frac{1-\delta(\rr)}{2} \rr^2
\]

\[
\rr (1 - \delta(\rr) - \zz(t)) \leq 2 \|D^{-1} \nabla L(\thetas) \|
\]

\[
\rr \leq 4 \|D^{-1} \nabla L(\thetas) \|
\]
From Assumptions 1, 2  follows that

\[
\| D (\opttheta - \theta^*) +  D^{-1} \{\nabla L (\opttheta) -  \nabla L (\thetas) \}   \|
\leq  \rombt
\]

\[
\| D (\opttheta - \theta^*) -  D^{-1} \nabla L (\thetas) \|
\leq  \rombt
\]
Not that for the coordinates transform $S$  there is an invariant:

\[
\left\|
\left(
\begin{matrix}
\breve{D} & 0 \\
0 & D_{\vartheta}
\end{matrix} \right)
\left( \begin{matrix}
u  - u^* \\
\vartheta - \vartheta^s
\end{matrix} \right)
\left(
\begin{matrix}
\breve{D}^{-1} & 0 \\
0 & D_{\vartheta}^{-1}
\end{matrix} \right)
\left( \begin{matrix}
\scorer L(u,  \vartheta)  - \scorer L(u^*,  \vartheta^*)\\
\nabla_\vartheta L(u,  \vartheta)  + \nabla_\vartheta L(u^*,  \vartheta^*)
\end{matrix} \right)
\right\|
\]

\[
= \|  D(\theta - \thetas) +  D^{-1} \{\nabla L (\theta) - \nabla L (\thetas) \} \|
\]

\[
\left\|
\left(
\begin{matrix}
\breve{D} & 0 \\
0 & D_{\vartheta}
\end{matrix} \right)
\left( \begin{matrix}
u  \\
\vartheta
\end{matrix} \right)
\right \| ^2  =  \theta ^T S^T [(S^{-1})^T D^2 (S^{-1})] S \theta  = \| D  \theta  \|^2
\]

\[
\left \| \left(
\begin{matrix}
\breve{D}^{-1} & 0 \\
0 & D_{\vartheta}^{-1}
\end{matrix} \right)
\left( \begin{matrix}
\scorer \\
\nabla_\vartheta
\end{matrix} \right)
\right\|^2 =
\nabla^T S^{-1}   [(S^{-1})^T D^2 (S^{-1})] ^{-1}  (S^{-1})^T \nabla =  \nabla^T D^{-2} \nabla
\]

\[
\left( \begin{matrix}
\scorer \\
\nabla_\vartheta
\end{matrix} \right)^T \left( \begin{matrix}
u  \\
\vartheta
\end{matrix} \right)  = \nabla^T S^{-1} S \theta = \nabla^T \theta
\]
Subsequently basing on this invariant

\[
\| \breve{D} (\widehat{u} - u^*) -  \breve{D}^{-1} \scorer L (\thetas) \|
\]

\[
\leq \| D (\opttheta - \theta^*) -  D^{-1} \nabla L (\thetas) \| \leq \rombt
\]

§.§ Independent models

Consider independent models (models with independent observations) and obtain a simpler variant of the Assumption 2 for this case. Involve three basic Lemmas for that.

Let $X_{1}\ldots X_{n}$ be independent real-valued random variables.
Assume that exist positive numbers $\mathbf{v}$ and $R$ such that

\[
\mathbf{v^{2}=}\sum_{i=1}^{n}\E X_{i}^{2}
\]

and for all integers $q\geq3$

\[
\sum_{i=1}^{n}\E\left[X_{i}\right]_{+}^{q}\leq\frac{q!}{2}\mathbf{v^{2}}R^{q-2}
\]

Then for all $\lambda\in(1,1/R)$

\[
\log\E e^{\lambda\sum_{i}(X_{i}- \E X_{i})}\leq\frac{\mathbf{v}^{2}\lambda^{2}}{2(1-R\lambda)}
\]

Let $\Omega$ be a finite pseudometric space and let $f(\theta)$ ($\theta\in \Omega$) be a collection of random variables such that for some constants $a,\mathbf{v},R>0$,
for all $\theta_1,\theta_2 \in \Omega$ and all $0 < \lambda < (Rd(\theta_1,\theta_2))^{-1}$

\[
\log\E \exp \{\lambda(f(\theta_1) -  f(\theta_2)) \}\le a\lambda d(\theta_{1}, \theta_{2})+\frac{\mathbf{v^{2}}\lambda^{2}d^{2}(\theta_{1}, \theta_{2})}{2(1 - R\lambda d(\theta_{1}, \theta_{2}))}
\]

Then for any $\theta_{0}\in \Omega$,

\[
\E[\sup_{\theta} f(\theta) - f(\theta_{0}) ]\leq3a\rr+12\mathbf{v}\int_{0}^{\rr/2}\sqrt{\log N(\varepsilon,\Omega)}d\varepsilon+12R\int_{0}^{\rr/2}\log N(\varepsilon,\Omega)d\varepsilon
\]

where $\rr=\sup_{\theta \in \Omega}d(\theta, \theta_{0})$ and $N(\varepsilon, \Omega)$ is covering number.

Consider independent random
variables $X_{1}\ldots X_{n}$ and let
be countable set of functions that satisfy conditions $\E f(X_{i})=0$
and $\|f\|_{\infty}\leq R$. Define

\[
Z=\sup_{f\in\mathcal{F}} \sum_{i=1}^{n}f(X_{i})
\]

Let $\textbf{v}^{2}\geq  \sup_{f \in \mathcal{F}}  \sum_{i=1}^{n}\E f^{2}(X_{i})$
then with probability $\hiprobt$

\[
Z<\E Z + \sqrt{2t(\textbf{v}^{2}+2R\E Z)}   +\frac{t R}{3}
\]

Apply three previous Lemmas in order to simplify  Assumption 2 for  independent models. Likelihood  of an independent model is a sum of independent functions:

\[
(L-\E L)(\theta) = \zeta(\theta) =  \sum_{i=1}^n \zeta_i(\theta)
\]

Note that $\zeta_i $ depends from the implicit i-th element from the dataset, such that $\zeta_i(\theta) = \zeta_i(\theta, Y_i)$.
Let $D^2$ be the Fisher matrix defined above and  $\forall \theta \in \Omega(\rr)$

\[
\sup_{\| u \| = 1} \sum_{i=1}^n \E (u^T D^{-1} \nabla^2 \zeta_i (\theta) D^{-1} u)^2 \leq \textbf{v}^2
\]

\[
\| D^{-1} \nabla^2 \zeta_i (\theta) D^{-1} \| \leq R
\]

Then Assumption 2 fulfills inside $\Omega(\rr)$ with probability $\hiprobt$ and

\[
\zz(t) \leq  \left( E +  \sqrt{2t(\textbf{v}^{2}+2R E)}   +\frac{tR}{3}  \right)
\]

\[
E = \frac{ 12 \mathbf{v} }{\rr } \int_{0}^{\infty}\sqrt{2 \log N(\varepsilon, \Omega(\rr))}d\varepsilon
+ \frac{12R}{\rr }   \int_{0}^{\infty}2 \log N(\varepsilon, \Omega(\rr) )d\varepsilon
\]

Set a random process for each $i$:

\[
X_i(\gamma, \theta)  = \frac{1}{\rr} \gamma^T \{ \nabla\zeta_i(\theta) -  \nabla\zeta_i(\thetas) \}
\]

Such that

\[
\sup_{\|D \gamma\| \leq \rr} \sum_i X_i(\gamma, \theta) = \| D^{-1} \{ \nabla\zeta(\theta) -  \nabla\zeta(\thetas) \} \|
\]

$\forall \, \text{fixed} \, (\gamma, \theta) \in \Omega(\rr, 0) \times \Omega(\rr, \thetas)$ and $\| u \| = 1$:

\[
\sup_u \E \sum_i (\nabla_{\theta} X_i(\gamma, \theta)^T D^{-1} u)^2
= \sup_u  \E \sum_i \frac{1}{\rr} ( \gamma \nabla^2 \zeta (\theta)^T D^{-1} u)^2
\]

\[
\leq \sup_u \E \sum_i ( u^T D^{-1} \nabla^2 \zeta (\theta)^T D^{-1} u)^2
\leq \textbf{v}^2
\]

\[
\sup_u \E \sum_i (\nabla_{\gamma} X_i(\gamma, \theta)^T D^{-1} u)^2 \leq \textbf{v}^2
\]

$\forall i \in {1, \ldots ,n} :$

\[
\| D^{-1} \nabla X_i(\gamma, \theta) \| \leq R
\]

Apply Lemma <ref> for the sum of random variables $X (\gamma, \theta) =  \sum_i X_i(\gamma, \theta)$ when $(\gamma, \theta)$ are fixed.

logexpλ( X(γ_1, θ_1) - X(γ_2, θ_2) )

= logexpλ( (γ_1 - γ_2)^T ∇_γ X(γ, θ) ) + logexpλ( (θ_1 - θ_2)^T  ∇_θ X(γ, θ) )

≤sup_u logexpλ(  D(γ_1 - γ_2)  u^T D^-1 ∇_γ X(γ, θ) )

+ sup_u logexpλ(  D (θ_1 - θ_2)  u^T D^-1  ∇_θ X(γ, θ) )

≤𝐯^2λ^2  D (γ_2 - γ_1) ^2  /2(1-Rλ D (γ_2 - γ_1)   ) + 𝐯^2λ^2  D(θ_2 - θ_1) ^2 /2(1-Rλ D(θ_2 - θ_1) )

≤𝐯^2λ^2 d^2_12 /2(1-Rλd_12)

\[
d^2_{12} = \| D(\theta_2 - \theta_1) \|^2  + \| D( \gamma_2 - \gamma_1) \|^2
\]

\[
\Upsilon = \Omega(\rr) \times \Omega(\rr)
\]

such that $\log N(\varepsilon, \Upsilon) = 2 \log N(\varepsilon, \Omega (\rr))$. Then with Lemma <ref> we obtain

\[
\E \sup_{\gamma, \theta} X(\gamma, \theta) \leq  12 \mathbf{v}  \int_{0}^{\infty}\sqrt{\log N(\varepsilon, \Upsilon)}d\varepsilon
+ 12R  \int_{0}^{\infty}\log N(\varepsilon, \Upsilon)d\varepsilon
\]

Applying Lemma <ref> to the random variable $Z = \sup_{\gamma, \theta} X(\gamma, \theta)$ completes the proof.

§.§ Covering numbers and entropy

Below one can read a short excerpt about an entropy of  ball and ellipsoid. The general formula for the covering number $N$ of a  convex set $\Omega$ in $\mathbb{R}^{p}$ with an arbitrary distance $d(\theta_1, \theta_2)$ is

\[
\]
where $B_1$ is a unit ball.

Ball entropy: Let $\Omega = B_\rr$ and $d(\theta_1, \theta_2) =  \| \theta_1 - \theta_2 \| $ then

\[
N(\varepsilon,  B_1) \leq \left(1+\frac{2}{\varepsilon}\right)^{p}
\]

and since $N( \varepsilon \rr, B_\rr) = N( \varepsilon, B_1)$

\[
\int_0^{\infty} \sqrt{\log N( \varepsilon, B_\rr)} d \varepsilon =
\rr  \int_0^{1} \sqrt{\log N( \varepsilon,  B_1 )} d \varepsilon
\]

\[
\leq \rr \sqrt{p} \int_{0}^1 \sqrt{ \log(3/\varepsilon) } d\varepsilon \leq  1.42 \rr  \sqrt{p}
\]

\[
\int_0^{\infty} \log N ( \varepsilon, B_\rr) d \varepsilon \leq
2.1 \rr  p
\]

Ellipsoid entropy: Let $\Omega = \mathcal{E}_\rr(D)$ and $d(\theta_1, \theta_2) =  \| D (\theta_1 - \theta_2) \| $. The entropy in this case is rather complicate in calculation. So we provide here only the the final statement from V. Spokoiny's  lecture notes [9].

\[
\int_0^{\infty} \sqrt{\log N( \varepsilon,  \mathcal{E}_\rr(D))} d \varepsilon
\lesssim \rr \sqrt{\alpha - 1} \sqrt{\sum_i \frac{ \log^\alpha (\lambda_i^2(D)) }{\lambda_i^{2}(D)}  }
\]

\[
\int_0^{\infty} \log N( \varepsilon, \mathcal{E}_\rr(D)) d \varepsilon
\lesssim \rr  \sum_i \frac{ 1 }{\lambda_i(D)}
\]

§ BARYCENTERS MODEL

We are going to show that Assumptions 1,2,3 are fulfilled for the barycenters model defined in Section <ref>. Also we need to estimate $\rombt$.  Remind that we deal with Likelihood function $L(\theta) = L(\theta, \{\theta_i\}_{i=1}^n)$ where  implicit random vectors $\{\theta_i\}_{i=1}^n$ is a dataset of Fourier coefficients corresponded to the  random measures $\{ \mu_i\}_{i=1}^n$.

Assumption 1:

\[
\| D^{-1} \{ \nabla^2 \E L(\theta) - \nabla^2 \E L(\thetas) \} D^{-1} \|
\leq  \| D^{-1} \{ \nabla^3 \E L(\theta) D^{-1} \} D^{-1} \| \rr
\]

Let $q(\theta_i)$ be the  distribution of each $\theta_i$  then

\[
\nabla^3 \E_i L(\theta - \theta_i) = \sum_{i=1}^n \int \nabla^3 l(\theta - \theta_i) q(\theta_i) d\theta_i = - \sum_{i=1}^n \int \nabla^2 l(\theta - \theta_i) \times \nabla q(\theta_i) d\theta_i
\]

\[
\| D^{-1} \{ \nabla^3 \E L(\theta) D^{-1} \} D^{-1} \| \leq \int \| D^{-1}  \nabla^2  L(\theta - \theta_x)  D^{-1} \|  \| D^{-1} \nabla p(\theta_x) \| d\theta_x
\]

and from the consequence of Theorem <ref> one gets

\[
\|  D^{-1} \nabla^2 l(\theta - \theta_i) D^{-1} \| \leq \frac{1}{ \varepsilon \lambda_{\min} \big(D K \circ G D \big) }
\]

\[
\| D^{-1} \{ \nabla^2 \E L(\theta) - \nabla^2 \E L(\thetas) \} D^{-1} \| \leq  \frac{\rr n }{ \varepsilon \lambda_{\min} \big(D K \circ G D \big)} \int  \| D^{-1} \nabla q(\theta) \| d\theta
\]

\[
\delta(\rr) =  \frac{\rr n }{ \varepsilon \lambda_{\min} \big(D K \circ G D \big)} \int  \| D^{-1} \nabla q(\theta) \| d\theta
\]

Assumption 2: From Theorem <ref> follows that if:

\[
\frac{  n }{ \varepsilon^2\lambda_{\min}^2 \big(D K \circ G D \big)} \leq \textbf{v}^2
\]

\[
\frac{ 1 }{\varepsilon \lambda_{\min} \big(D K \circ G D \big)} \leq R
\]

\[
\zz(t) \leq  E +  \sqrt{2t(\textbf{v}^{2}+2R E)}   +\frac{tR}{3}
\]

\[
E = 12 \mathbf{v} \sqrt{ 2 p_D} + 24 R \,   p_D
\]

and $p_D$ is ellipsoid entropy with matrix $D$

\[
p_D = \sqrt{\sum_i \frac{ \log^2 (\lambda_i^2(D)) }{\lambda_i^{2}(D)}  }
\]

Assumption 3: Each model component $l(\theta - \theta_i)$ is convex since

\[
l(\lambda \theta_1 + (1 - \lambda) \theta_2 - \theta_i) =  l(\lambda (\theta_1 - \theta_i) + (1 - \lambda) (\theta_2 - \theta_i))
\]

\[
= \max_{\eta \in W_{2,1}} \langle \eta, \lambda (\theta_1 - \theta_i) + (1 - \lambda) (\theta_2 - \theta_i)) \rangle
\]

\[
\leq \max_{\eta \in W_{2,1}} \langle \eta, \lambda (\theta_1 - \theta_i) \rangle + \max_{\eta \in W_{2,1}} \langle \eta, (1 - \lambda) (\theta_1 - \theta_i) \rangle
\]

\[
= \lambda l(\theta_1 - \theta_i) + (1 - \lambda) l(\theta_2 - \theta_i)
\]

Note that $l^2$ is also convex as a composition of convex functions and the complete model $L$ is convex $(\nabla^2 L > 0)$ as a positive aggregation of convex functions.  Combination of these assumptions is used in the proof of Theorem <ref> which gives the deviation bound of parameter $\widehat{\theta}$.

§ SUPPORT FUNCTIONS

Bounds for the first and second derivatives of the Likelihood of barycenters model involves additional theory from Convex analysis.
Legendre–Fenchel transform of a function $f:X \to \overline{\R}$ or the convex conjugate function calls

\[
f^*(y) = \sup_{x \in X} ( \langle x, y \rangle - f(x) )
\]

Support function for a convex body $E$ is

\[
s(\theta) = \sup_{\eta \in E}  \theta^T \eta
\]

Note that for indicator function $\delta_E(\eta)$ of a convex set $E$ the conjugate function is support function of $E$

\[
\delta^*_E(\theta) = s(\theta)
\]

[$\oplus$]  Let $f_1, f_2:E \to \overline{\R}$ are convex functions. The infimal convolution of them is

\[
(f_1 \oplus f_2)(x) = \inf_{x_1 + x_2 = x} (f_1(x_1) + f_2(x_2))
\]

[2] Let $f_1, f_2:E \to \overline{\R}$ are convex lower-semi-continuous functions. Then

\[
(f_1 \oplus f_2)^* = f_1^* + f_2^*
\]

\[
(f_1 + f_2)^* = \overline{f_1^* \oplus f_2^*}
\]

The support function of intersection $E = E_1 \cap E_2$ is infimal convolution of support functions for $E_1$ and $E_2$

\[
s(\theta) = \inf_{\theta_1 + \theta_2 = \theta} (s_1(\theta_1) + s_2(\theta_2))
\]

\[
s_1(\theta) = \sup_{\eta \in E_1}  \theta^T \eta,  \quad
s_2(\theta) = \sup_{\eta \in E_2}  \theta^T \eta
\]

\[
\delta_{E_1 \cap E_2}(\eta) = \delta_{E_1}(\eta) + \delta_{E_2}(\eta),
\]

\[
(\delta_{E_1} + \delta_{E_2})^* = \overline{ \delta^*_{E_1} \oplus \delta^*_{E_2} }
\]

With additional property

\[
\text{intdom} \, \delta_{E_1} \cap \text{dom} \, \delta_{E_2} = \text{int} E_1 \cap E_2 \neq  \emptyset
\]

one have

\[
(\delta_{E_1} + \delta_{E_2})^* = \delta^*_{E_1} \oplus \delta^*_{E_2}
\]

Let a support function $s(\theta)$ is differentiable, then its gradient lies on the border of corresponded convex set $E$

\[
\nabla s(\theta) =  \widehat{\eta}(\theta)  \in \partial E
\]

\[
\widehat{\eta}(\theta) = \argmax_{\eta \in E} \eta^T \theta
\]

It follows from the convexity of $E$ and linearity of optimization functional.

\[
\frac{\partial \widehat{\eta}(\theta) }{\partial \| \theta \| }  = 0 \Rightarrow \frac{\partial \widehat{\eta}(\theta) }{\partial  \theta  }^T \theta = 0
\]

\[
\nabla s(\theta) = \frac{\partial \widehat{\eta}(\theta) }{\partial  \theta  }^T \theta  + \widehat{\eta}(\theta) = \widehat{\eta}(\theta)
\]

Optimization related to support function.

[2]
Let $f_1, f_2:E \to \overline{\R}$ are convex continuous functions. Then  the subdifferential of their infimal convolution can be computed by following formula

\[
\partial (f_1 \oplus f_2)(x) = \bigcup_{x = x_1 + x_2} \partial f(x_1) \cap \partial f(x_2)
\]

If in addition $f_1, f_2$ are differentiable, then their infimal convolution is differentiable and $\exists x_1, x_2: x = x_1 + x_2$

\[
\nabla (f_1 \oplus f_2)(x) = \nabla f_1(x_1) = \nabla f_2(x_2)
\]

Let $f_1,\ldots, f_m:E \to \overline{\R}$ are convex and two times differentiable functions. There are following upper bound for the second derivative of the infimal convolution $\forall t: \sum_{i=1}^m t_i = 1$

\[
\partial \nabla^T (f_1 \oplus \ldots \oplus f_m)(x) \preceq \sum_{i=1}^m t_i^2 \nabla^2 f(x_i)
\]

where $\sum_{i=1}^m x_i  = x$.

Use notation $f = f_1 \oplus \ldots \oplus f_m$. Let

\[
f(y) = \sum_i f_i(y_i)
\]

According to Lemma <ref> if all the functions are differentiable then

\[
\nabla f(y) = \sum_i t_i \nabla f_i(y_i)
\]

From the definition $\oplus$ also follows that

\[
f(y + z) \le \sum_i f_i(y_i + t_i z)
\]

Make Tailor expansion for the left and right parts and account equality of the first derivatives.

\[
z^T \partial \nabla^T f(y + \theta z) z
\leq \sum_i  t_i^2 z^T \nabla^2 f_i(y_i + \theta_i z) z
\]

Since direction $z$ was chosen arbitrarily then dividing both parts of the previous equation by $\|z\|^2 \to 0$ we come to inequality

\[
\partial \nabla^T f(y)   \preceq   \sum_i  t_i^2   \nabla^2 f_i(y_i)
\]

One can find another provement of the similar Theorem in book [2] (Theorem 18.15).

Let $f_1,\ldots, f_m:E \to \overline{\R}$ are convex and two times differentiable functions. There are following upper bounds for infimal convolution $f = f_1 \oplus \ldots \oplus f_m$ derivatives $\forall \gamma$ and some matrix $A$

\[
\gamma^T \partial \nabla^T f(x) \gamma  \leq \max_i \gamma^T \nabla^2 f_i(x_i) \gamma \frac{  f(x_i)  }{  f(x) }
\]

\[
\gamma^T \partial \nabla^T f^2(x)\gamma \leq 2 (\gamma^T \nabla f(x))^2  +  2  \max_i \gamma^T \nabla^2 f_i(x_i) \gamma  f(x_i)
\]

Choosing appropriate $\{t_i\}$ in Lemma <ref> one get the required upper bounds. Set

\[
t_i = \frac{f(x_i)}{\sum_{j=1}^m f (x_j) }
\]

and since

\[
\sum_{j=1}^m f(x_j) = f(x)
\]

\[
\sum_i  t_i^2  \gamma^T \nabla^2 f_i(y_i) \gamma  \leq \max_i t_i \gamma^T \nabla^2 f_i(y_i) \gamma = \max_i \gamma^T \nabla^2 f_i(x_i) \gamma f(x_i)
\]

To prove the second formula apply this inequality in

\[
\partial \nabla^T f^2 = 2 \nabla f \nabla^T f + 2f \partial \nabla f
\]

Let $s_1,\ldots, s_m: E^* \to \overline{\R}$ are support functions of the bounded convex smooth sets $E_1, \ldots, E_m$. There is upper bound for the  derivatives of support function $s$ of intersection $E_1 \cap \ldots \cap E_m$, such that  $\forall i$

\[
\gamma^T \partial \nabla^T s(x) \gamma \leq \frac{ \max_i \gamma^T \partial \eta_i / \partial x_i  \gamma s_i (x_i)  }{  s ( x ) }
\]

\[
\gamma^T \partial \nabla^T s^2(x) \gamma \leq 2 (\gamma^T \eta_i)^2 +  2  \max_i \gamma^T \partial \eta_i / \partial x_i \gamma  s_i( x_i)
\]

It follows from Theorem <ref>  and Lemma <ref>.

§ WASSERSTEIN DISTANCE AS A SUPPORT FUNCTION

[W-dual]Consider two random variables $X$ and $Y\in\mathbb{R}^{p}$ with
densities $\varphi_{X}$ and $\varphi_{Y}$.  Define Wasserstein distance in dual form between them as

\[
W_1(\varphi_X, \varphi_Y)=\max_{\forall x: \|\nabla f(x)\|\leq 1}\{\E f(X)-\E f(Y)\}
\]

where $\forall x: \|\nabla f(x)\|\leq 1 $ means that function $f$ is $1$-Lipshits.

Note that if $\pi(x,y)$ is a joint distribution with marginals $\varphi_{X}$ and $\varphi_{Y}$ then this definition is equivalent to  the original definition

\[
W_1(\varphi_ X, \varphi_ Y)=\min_{\pi} \E \| X - Y \|
\]

which follows from Kantorovich-Rubinstein duality [6].  Involve a normalized Fourier basis $\left\{ \psi_{k}(x)\right\} _{k\in\mathbb{N}^{p}}$
with scalar product  Gram function $G(x)$.
[$W-dual-regularised$]Consider two random variables $X$ and $Y\in\mathbb{R}^{p}$ with
densities $\varphi_{X}$ and $\varphi_{Y}$. Define a penalized Wasserstein distance
between them in dual form as

\[
\widetilde{W}(\varphi_ X, \varphi_ Y)=\max_{\forall x: \|\nabla f(x)\|\leq 1} \left\{\E f(X)-\E f(Y) -   \varepsilon \int  \| \nabla f (x) \|^2  G(x) dx   \right\}
\]

The regulariser term in this definition allows to bound the second derivative of the distance which will be shown below.
Wasserstein distance is a support function (ref. Def(s)) in Fourier basis.  In which connection

\[
\]

\[
\eta_{k}(f)=\langle f,\psi_{k}\rangle_{G}=\int f(x)\psi_{k}(x)G(x)dx
\]

Now we can rewrite the expectation difference as

\[
Ef(X)-Ef(Y)=\langle f,\frac{\varphi_{X}}{G}\rangle_{G}-\langle f,\frac{\varphi_{Y}}{G}\rangle_{G}=\langle\eta(f),\theta(\varphi_{X})\rangle-\langle\eta(f),\theta(\varphi_{Y})\rangle
\]

\[
\theta_{k}(\varphi)=\int\varphi(x)\psi_{k}(x)dx
\]

Define positive symmetric matrices

\[
K_x =  \left(\begin{matrix}
\nabla^T \psi_1(x) \\
\ldots \\
\nabla^T \psi_k(x) \\
\ldots
\end{matrix} \right)
\left(\begin{matrix}
\nabla \psi_1(x) & \ldots & \nabla \psi_k(x) & \ldots
\end{matrix} \right) = (\nabla^T \psi(x)) (\nabla \psi^T(x))
\]

\[
K \circ G = \int K_x G(x) dx
\]

Each $K_{x}$ is positive, since $\eta^{T}K_{x}\eta=\|\nabla f(x)\|^{2}$.
Condition $\forall x: \|\nabla f(x)\|\leq 1$ is equivalent
in Fourier basis to

\[
\eta \in  \bigcap \mathcal{E}_x =\left\{ \eta : \forall x: \left( \sum_{k}\eta_{k}\nabla\psi_{k}(x) \right)^{2}=\eta^{T} K_x \eta \leq 1 \right\}
\]

An important remark is that

\[
\bigcap \mathcal{E}_x  \subset \left\{ \eta : \eta^{T} (K \circ G) \eta \leq 1 \right\}
\]

Finally we have come to the Wasserstein distance in Fourier basis.

Let random vectors $X$ and $Y$ have
densities $\varphi_{X}$ and $\varphi_{Y}$ with Fourier coefficients $\theta_X$ and $\theta_Y$, then the Wasserstein distance  is the support function of the convex set $\bigcap \mathcal{E}_x$ defined above, i.e.

\[
W_1( \varphi_X,  \varphi_Y)= \max_{\eta \in \bigcap \mathcal{E}_x }\langle\eta,\theta_X -\theta_Y \rangle
\]

As for regularised case

\[
\widetilde{W_1} (\varphi_X, \varphi_Y)= \max_{\eta \in \bigcap \mathcal{E}_x }\langle\eta,\theta_X - \theta_Y\rangle - \varepsilon \eta^T  K \circ G \eta
\]

Remind that barycenters Likelihood consists of independent components $l_i(\theta - \theta_i)$ with a random vectors $\theta_i \in \R^{\infty}$ and parameter $\theta \in \R^{\infty}$.

\[
l(\theta - \theta_i) = \max_{\eta \in \bigcap \mathcal{E}_x }\langle\eta,\theta - \theta_i \rangle
- \varepsilon \eta^T  K \circ G \eta
\]
Note that by definition  the dual function of $l$ is

\[
l^*(\eta)  = \delta_{\bigcap \mathcal{E}_x}(\eta) + \varepsilon \eta^T  K \circ G \eta
\]

Consequently from Lemma <ref> follows that

\[
l(\theta - \theta_i) = \delta^*_{\bigcap \mathcal{E}_x}(\theta - \theta_i) \oplus (\varepsilon \eta^T  K \circ G \eta)^*(\theta - \theta_i)
\]

\begin{equation}
\label{l_conv}
=\max_{\eta \in \bigcap \mathcal{E}_x }\langle\eta,\theta - \theta_i \rangle
\oplus \frac{1}{\varepsilon} (\theta - \theta_i)^T  (K \circ G)^{-1} (\theta - \theta_i)
\end{equation}

Apply Theorem <ref>  taking into account $\bigcap \mathcal{E}_x \subset\left\{ \eta : \eta^{T} (K \circ G) \eta \leq 1 \right\}$, one gets the following bounds on the  derivatives of the function $l$.
The gradient upper bounds for functions $l$ and $l^2$ are

\[
\| D^{-1} \nabla l \| \leq \frac{1}{\lambda^{1/2}_{\min} \big( D  K\circ G D \big)}
\]

\[
\| D^{-1} \nabla l^2 (\theta - \theta_i) \| \leq \frac{ 2 \| (K\circ G)^{-1/2} (\theta - \theta_i) \| }{\lambda^{1/2}_{\min} \big( D  K\circ G D\big)}
\]

\[
\etas(\theta) = \argmax_{\eta \in \bigcap \mathcal{E}_x}  \eta^T \theta
\]

Use the equation (<ref>). By the consequence of Lemma <ref> and Lemma <ref> $\exists \theta_0$:

\[
\nabla l (\theta - \theta_i)  = \eta^*(\theta_0)
\]

Since $\| (K\circ G)^{1/2} \etas \| \leq 1$

\[
\| D^{-1} \nabla l \| = \| D^{-1} \etas \|  = \| D^{-1} (K\circ G)^{-1/2} (K\circ G)^{1/2} \etas \|  \leq \| D^{-1} (K\circ G)^{-1/2} \|
\]

and from $ \nabla l^2 = 2 l \nabla l$ one gets

\[
\| D^{-1} \nabla l^2 (\theta - \theta_i) \| \leq 2 l(\theta - \theta_i) \| D^{-1} \nabla l \| \leq 2 \| (K\circ G)^{-1/2} (\theta - \theta_i) \| \| D^{-1} \nabla l \|
\]

The second derivative  upper bounds for functions $l$ and $l^2$ are

\[
\|  D^{-1} \partial \nabla^T l(\theta - \theta_i)  D^{-1} \| \leq
\frac{1}{\min_x \lambda_{\min}(D  K_x D) \| (K\circ G)^{-1/2} (\theta - \theta_i) \|}
\]

\[
\|  D^{-1} \partial \nabla^T l(\theta - \theta_i)  D^{-1} \| \leq
\frac{1}{ \varepsilon \lambda_{\min}(D  K \circ G D)  }
\]

\[
\|  D^{-1} \partial \nabla^T l^2  D^{-1} \| \leq
\frac{2}{\min_x \lambda_{\min}(D K_x D)}
\]

Matrix $K_x$ may be singular which makes the first bound non-informative. The second bound comes from the regulariser  $\varepsilon \eta^T  K \circ G \eta $ and has big coefficient $(1 / \varepsilon)$.  It is a weak part of the current theory and requires an improvement or an example which shows that this bound it tight.

Consider support function with one ellipsoid.

\[
s_x (\theta) =  \max_{\eta^{T} K_x \eta\leq 1}\langle\eta,\theta\rangle=\|K_x^{-1/2}\theta\|
\]

Denote $\eta^{*}(\theta) = \argmax \langle\eta,\theta\rangle $, $\eta^{T} K_x \eta\leq 1$.

\[
\eta^{*}(\theta)=\frac{K_x^{-1}\theta}{\|K_x^{-1/2}\theta\|}
\]

\[
\frac{\partial\eta^{*}(\theta)}{\partial\theta}=\frac{K_x^{-1}\theta^{T}K_x^{-1}\theta-K_x^{-1}\theta\theta^{T}K_x^{-1}}{\left(\theta^{T}K_x^{-1}\theta\right)^{3/2}}
\]

For some vector $\| \gamma \| = 1$ and property $ \| a \|^2 \| b \|^2 \geq (a^Tb)^2 $

\[
\gamma^TK_x^{-1} \gamma \theta^{T}K_x^{-1}\theta - \gamma^T K_x^{-1}\theta\theta^{T}K_x^{-1} \gamma \leq \| K_x^{-1} \| \theta^{T}K_x^{-1}\theta
\]

\[
\left\| \frac{\partial\eta^{*}(\theta)}{\partial\theta} \right \| \leq \frac{\|K_x^{-1} \|}{\left(\theta^{T}K_x^{-1}\theta\right)^{1/2}}
\]

Apply Theorem <ref>

\[
\| D^{-1} \partial \nabla^T l(\theta - \theta_i) D^{-1} \| \leq \max_x \left \|D^{-1}\frac{\partial\eta_x^{*}(\thetas_x)}{\partial\theta} D^{-1} \right \| \frac{s_x(\thetas_x)}{s(\theta - \theta_i)}
\leq \frac{\max_x \|D^{-1} K_x^{-1} D^{-1}\|}{ \| (K\circ G)^{-1/2} (\theta - \theta_i) \| }
\]
The second bound for this norm follows directly from Lemma <ref> and equation (<ref>).
Now consider the squared Wasserstein distance ($l^2$) which has a better derivative bound. From Theorem <ref> one gets

\[
\|  D^{-1} \partial \nabla^T l^2 D^{-1}  \| \leq 2 \max_x \left \| D^{-1} \eta^{*}(\thetas_x) \eta^{*}(\thetas_x)^T D^{-1} +   D^{-1} \frac{\partial\eta^{*}(\thetas_x)}{\partial\theta} \|K_x^{-1/2}\thetas_x\| D^{-1} \right \|
\]

Note that

\[
\frac{\partial\eta^{*}(\theta)}{\partial\theta} \|K_x^{-1/2}\theta\| =
K_x^{-1}  - \frac{(K_x^{-1} \theta)(K_x^{-1} \theta)^T}{ \|K_x^{-1/2}\theta\|^2 }
\]

\[
\eta^{*}(\theta) \eta^{*}(\theta)^T +   \frac{\partial\eta^{*}(\theta)}{\partial\theta} \|K_x^{-1/2}\theta\| = K_x^{-1}
\]

\[
\|  D^{-1} \partial \nabla^T l^2 D^{-1}  \| \leq 2 \max_x \| D^{-1} K_x^{-1} D^{-1} \|
\]

Note that the Wasserstein distance also may be differentiated directly. In paper [8] one may find the lemma about directional derivative. For directions $h_1, h_2$ it holds

\[
d'_W(\mu_1, \mu_2)(h_1, h_2)  =  \max_{(u,v) \in \Phi(\mu_1, \mu_2)} - (\langle u,h_1 \rangle + \langle v,h_2 \rangle)
\]

\[
\Phi = \{ (u,v): \langle u, \mu_1 \rangle + \langle v, \mu_2 \rangle = d_W(\mu_1, \mu_2), \, \forall (x,y):
u(x) + v(y) \leq \| x - y \|  \}
\]

§ GAUSSIAN APPROXIMATION

The multivariate Hermite polynomial $H_k$
is defined by

\[
H_k(x) = (-1)^{|k|} e^{x^T \Sigma^{-1} x /2} \frac{\partial^{|k|}}{\partial^{k_1} \ldots \partial^{k_p}} e^{-x^T \Sigma^{-1} x/2}
\]

Consider a Gaussian vector $Z \sim \ND(0, \Sigma)$ and two functions $h$ and $f_h$ such that

\[
f_h(x) = - \int_0^1 \E \overline{h}(Z(x,t))dt
\]

\[
\overline{h}(Z(x,t)) = h(\sqrt{t}x + \sqrt{1-t} Z) - \E h(Z)
\]

Then $f_h$ is a solution of the Stein equation

\[
\overline{h}(x) = (\tr\{\nabla^2 \Sigma \} - x^T \nabla) f_h(x)
\]

\[
\frac{\partial^{|k|}}{\partial^{k_1} \ldots \partial^{k_p}} f_h(x) =
- \int_0^1 \frac{1}{2} \frac{t^{ \frac{|k|}{2}  - 1}}{(1 - t)^{\frac{|k|}{2}  - 1}} \E H_k(Z) \overline{h}(Z(x,t)) dt
\]

\[
\nabla^2 f_h(x) - \nabla^2 f_h(y) =
- \int_0^1  \frac{1}{2(1 - t)} \E H_2(Z) \{ h(Z(x,t)) - h(Z(y,t)) \}  dt
\]

\[
H_2(Z) = (\Sigma^{-1}Z) (\Sigma^{-1}Z)^T  - \Sigma^{-1}
\]

\[
=  \Sigma^{-1/2} \{ (\Sigma^{-1/2}Z) (\Sigma^{-1/2}Z)^T  - I\} \Sigma^{-1/2}
\]

\[
X = \sum_{i=1}^n X_i
\]

\[
\E \overline{h}(X) =   \E \tr\{\nabla^2 \Sigma \} f_h(X) - \E \sum_{i=1}^n X_i^T \nabla  f_h(X)
\]

\[
= \sum_{i=1}^n \E X_i^T \left \{ \nabla^2 f_h(X' + \theta (X_i - X'_i) ) - \nabla^2 f_h(X') \right \} (X_i - X'_i)
\]

\[
= \sum_{i=1}^n \E (\Sigma^{-1/2}X_i)^T \Sigma^{1/2} \left \{ \nabla^2 f_h(X' + \theta (X_i - X'_i) ) - \nabla^2 f_h(X') \right \} \Sigma^{1/2} (\Sigma^{-1/2} (X_i - X'_i))
\]

For a unit vector $\| \gamma \| = 1$ and conditional expectation $\E_{-i} =  \E (\cdot | X_i, X'_i)$

\[
\gamma^T \E_{-i} \Sigma^{1/2} \left \{ \nabla^2 f_h(X' + \theta (X_i - X'_i) ) - \nabla^2 f_h(X') \right \} \Sigma^{1/2} \gamma
\]

\[
=  \int_0^1  \frac{1}{2(1 - t)} \E_{-i} \{ ((\Sigma^{-1/2}Z)^T \gamma)^2  - 1 \} \{ h(Z(X' + \theta (X_i - X'_i),t)) - h(Z(X',t)) \}  dt
\]

\[
\leq \int_0^{1-\alpha} \frac{t^{1/2}}{2(1 - t)} A dt
+  \int_{1-\alpha}^1 \frac{1}{(1 - t)^{1/2}} B dt
\]

\[
\leq - \frac{A}{2} \log(\alpha) + 2B \sqrt{\alpha}
\]

\[
\leq A \left(1 + \log \left(\frac{2B}{A} \right) \right)
\]

\[
A = \| X_i - X'_i \| \, \E_{-i} | ((\Sigma^{-1/2}Z)^T \gamma)^2  - 1 | \, \|  \nabla h( \sqrt{t} (X' + \theta_1 (X_i - X'_i) + \sqrt{1-t} Z ) \|
\]

\[
B = \E_{-i} | ((\Sigma^{-1/2}Z)^T \gamma)^2  - 1 | \, \| Z \| \, \|  \nabla h( \sqrt{t} (X' + \theta (X_i - X'_i) + \theta_2 \sqrt{1-t} Z ) \|
\]

Consider a sequence of independent zero-mean random vectors $X = \sum_{i=1}^n X_i$ in $\R^{p}$ with a covariance matrix

\[
\E X X^T  =  \Sigma
\]

Then the  Wasserstein distance between $X$ and Gaussian vector $Z \in \ND(0, \Sigma)$ has following upper bound

\[
d_W(X, Z) \leq \sqrt{2} \mu_3 \left (1 + \log(2 \sqrt{ \tr \{ \Sigma \} } \mu_2) - \log (\mu_3) \right)
\]

\[
\mu_3 =  \sum_{i=1}^n  \E \| \Sigma^{-1/2} (X_i - X'_i) \| \| \Sigma^{-1/2} X_i \| \| X_i - X'_i \|
\]

\[
\mu_2 =  \sum_{i=1}^n  \E \| \Sigma^{-1/2} (X_i - X'_i) \| \| \Sigma^{-1/2} X_i  \|
\]

In i.i.d case with $\Sigma = I_p$

\[
d_W(X, Z)  = O \left( \frac{p^{3/2} \log(n)}{\sqrt{n}} \right)
\]

These is the same theorem with a different provement in paper [3].

Consider a sequence of independent zero-mean random vectors $X = \sum_{i=1}^n X_i$ in $\R^{p}$ with  a covariance matrix

\[
\E X X^T  =  \Sigma
\]

Let a function $\varphi: \R^{p} \to \R_{+}$ be sub-additive:

\[
\varphi(x + y) \leq \varphi(x) + \varphi(y)
\]

and  with Gaussian vector $Z \in \ND(0, \Sigma)$ fulfills the anti-concentration property, such that

\[
\P (\varphi(Z) > x) - \P (\varphi(Z) > x + \Delta) \leq C_A \Delta
\]

Then the measure difference between $X$ and Gaussian vector $Z$ has following upper bound $\forall x$

\[
\left | \P (\varphi(X) > x) -  \P (\varphi(Z) > x) \right | \leq 22 C_A \mu_3 \log \left( \frac{4 \mu_2}{C_A \mu_3} \right) \log \left( \frac{\sqrt{2 \E \varphi^2(Z)}  \mu_2}{ 20 C_A \mu^2_3} \right)
\]

\[
\mu_3 =  \sum_{i=1}^n  \E \| \Sigma^{-1/2} (X_i - X'_i) \| \| \Sigma^{-1/2} X_i  \| \varphi( X_i - X'_i )
\]

\[
\mu_2 =  \sum_{i=1}^n  \E \| \Sigma^{-1/2} (X_i - X'_i) \| \| \Sigma^{-1/2} X_i  \|
\]

Define a smooth indicator function

\[
g_{x,\Delta}(t) = \begin{cases}
0, & t < x \\
(t-x)/\Delta, & t \in [x, x + \Delta] \\
1, & t > x + \Delta
\end{cases}
\]

Set $h = g_{x, \Delta} \circ \varphi$.
Denote the required bound by $\delta$:

| (φ(X) > x) -  (φ(Z) > x) |

≤max_Δ | g_x,Δ ∘φ(X) - g_x,Δ ∘φ(Z) |  ≤δ
Note that from sub-additive property of the function $\varphi$ follows

\[
g_{x, \Delta} \big ( \varphi(X + dX) \big) \leq g_{x, \Delta}  \big( \varphi(X) + \varphi( d X) \big)
\]

\[
g'_{x, \Delta}(t) = \frac{1}{\Delta} \Ind [ x < t < x + \Delta ]
\]

\[
\E g'_{x, \Delta}(\varphi(Z)) = \frac{1}{\Delta} \big( \P (\varphi(Z) > x) - \P (\varphi(Z) > x + \Delta) \big) \leq C_A
\]

\[
\E g'_{x, \Delta}(\varphi(Z(X,t))) \leq \frac{1}{\Delta} \big( \P (\varphi(Z) > x) - \P (\varphi(Z) > x + \Delta) \big) + \frac{2\delta}{\Delta} \leq C_A + \frac{2\delta}{\Delta}
\]

_-i h(Z(X' + θ(X_i - X'_i),t))  - _-i h(Z(X',t))

≤_-i g_x, Δ  ( φ(Z(X',t)) + φ( X_i - X'_i ) ) - _-i g_x, Δ  (  φ(Z(X',t) )

≤_-i g'_x, Δ ( φ(Z(X',t)) + θφ( X_i - X'_i ) )  φ( X_i - X'_i )

≤( C_A + 2δ/Δ ) φ( X_i - X'_i )

\[
\E h(Z(X',t)) - \E h(Z(X' + \theta (X_i - X'_i),t)) \leq \left ( C_A + \frac{2\delta}{\Delta} \right ) \varphi( X_i - X'_i )
\]

Apply this inequality denoting $\varepsilon^2 = (\Sigma^{-1/2}Z)^T \gamma)^2 \sim \ND^2(0, 1)$

_-i { ((Σ^-1/2Z)^T γ)^2  - 1 } { h(Z(X' + θ(X_i - X'_i),t)) - h(Z(X',t)) }

≤_-i  ε^2  { g_x, Δ   [ φ(Z(X',t)) + φ( X_i - X'_i )   ] - g_x, Δ   [ φ(Z(X',t))  ] }

+ _-i h(Z(X',t)) -   _-i h(Z(X' + θ(X_i - X'_i),t))

≤| τ- 1 | ( C_A + 2δ/Δ ) φ( X_i - X'_i ) +  [ε^2 > τ] ε^2

Let a random variable $\varepsilon$ has a tail bound $\forall \xx \geq \xx_0$

\[
\P(\varepsilon > h(\xx)) \leq e^{-\xx}
\]

Then for a function $g: \R_+ \to \R_+$ with derivative $g': \R_+ \to \R_+$

\[
\E \Ind[\varepsilon > h(\xx_0)] g(\varepsilon) \leq g(h(\xx_0)) e^{-\xx_0} + \int_{\xx_0}^{\infty} e^{-\xx} g'(h(\xx)) h'(\xx) d\xx
\]

In particular

\[
\E \Ind[\varepsilon > h(\xx_0)] \varepsilon \leq h(\xx_0) e^{-\xx_0} + \int_{\xx_0}^{\infty} e^{-\xx} h'(\xx) d\xx
\]

\[
\E \Ind[\varepsilon > h(\xx_0)] \varepsilon^r \leq h(\xx_0)^r e^{-\xx_0} + r \int_{\xx_0}^{\infty} e^{-\xx} h(\xx)^{r-1} h'(\xx) d\xx
\]

For $\varepsilon \sim \ND(0, 1)$ we have

\[
\P(\varepsilon > \sqrt{2\xx}) \leq e^{-\xx}
\]

and by means of the previous lemma we get

\[
\E \Ind[\varepsilon^2 > \tau] \varepsilon^2 = 2 \E \Ind[\varepsilon > \sqrt{\tau}] \varepsilon^2 \leq 2 (\tau + 2) e^{-\tau / 2}
\]

_-i { ((Σ^-1/2Z)^T γ)^2  - 1 } { h(Z(X' + θ(X_i - X'_i),t)) - h(Z(X',t)) }

≤|τ- 1| ( C_A + 2 δ/Δ ) φ( X_i - X'_i ) +   2 (τ+ 2) e^-τ/ 2

We need also another upper bound for this expectation when $t$ close to $1$.

_-i { ((Σ^-1/2Z)^T γ)^2  - 1 }   h(Z(X',t))

= _-i { ((Σ^-1/2Z)^T γ)^2  - 1 }  { h( √(t) X' +  √(1 - t) Z  ) - h( √(t) X') }

≤_-i | ((Σ^-1/2Z)^T γ)^2  - 1 |   | g'_x,Δ |   φ(√(1-t)Z)

≤1/Δ √(2 φ^2(Z))

Set $ \Delta = \delta / (2 C_A )$

\[
B =  \frac{2 C_A }{\delta} \sqrt{2 \E \varphi^2(Z)} \, \mu_2
\]

Set $\tau = 2 \log ( 4 \mu_2 / (C_A \mu_3) )$

A = 5 |τ- 1|  C_A  μ_3 +   2 (τ+ 2) e^-τ/ 2 μ_2

≤11 C_A μ_3 log( 4 μ_2/C_A μ_3 )

δ≤- A/2 log(α) + 2 B  √(α) + C_A Δ

≤2 A (1 + log(2B δ) - log(δ) - log(A) )

≤2 A (1 + log(2B δ)  - 2 log(A) + loglog(2B δ) - loglog(A) )

≤22 C_A μ_3 log( 4 μ_2/C_A μ_3 ) log( √(2 φ^2(Z))  μ_2/ 20 C_A μ^2_3 )

In i.i.d case with $\Sigma = I_p$ and $\phi(x) = O ( \| x \| )$

\[
\left | \P (\varphi(X) > x) -  \P (\varphi(Z) > x) \right |  = O \left(  C_A \mu_3 \log^2(n) \right)
\]

Note that  lemma <ref>  improves  the classical Multivariate Berry–Esseen Theorem [4] for the case of sub-additive functions $\phi(x) = O ( \| x \| )$. Namely it answers the open question “Whether one can remove or replace the factor $p^{1/4}$ by a better one (eventually by 1)”.

In i.i.d case with $\Sigma = I_p$ and $\phi(x) = \| x \| $, which is rather common is statistical learning theory, one have $C_A = O(1/\sqrt{p})$ and

\[
\left | \P (\|X\| > x) -  \P (\|Z\| > x) \right |  = O \left( p \frac{  \log^2(n)}{\sqrt{n}} \right)
\]

The author thanks Prof. Roman Karasev, Prof. Vladimir Spokoiny and Prof. Dmitriy Dylov for discussion and contribution to this paper.

[1]
Martial Agueh and Guillaume Carlier.
Barycenters in the Wasserstein space.
SIAM Journal on Mathematical Analysis, 43(2):904–924, 2011.

[2]
Heinz H. Bauschke and Patrick L. Combettes.
Convex Analysis and Monotone Operator Theory in Hilbert Spaces.
Springer Publishing Company, Incorporated, 1st edition, 2011.

[3]
V. Bentkus.
A new method for approximations in probability and operator theories.
Lithuanian Mathematical Journal, 43(4):367–388, Oct 2003.

[4]
V. Bentkus.
On the dependence of the berry–esseen bound on dimension.
Journal of Statistical Planning and Inference, 2003.

[5]
Massart P. Boucheron S., Lugosi G.
Concentration inequalities: A nonasymptotic theory of independence.
Oxford University Press, 2013.

[6]
D.A. Edwards.
On the kantorovich–rubinstein theorem.
Expositiones Mathematicae, 29(4):387 – 398, 2011.

[7]
Nicolas Papadakis Jérémie Bigot, Elsa Cazelles.
Penalized barycenters in the wasserstein space.
arXiv:1606.01025, 2016.

[8]
Axel Munk Max Sommerfeld.
Inference for empirical wasserstein distances on finite spaces.
arXiv:1610.03287v2, 2016.

[9]
V. Spokoiny.
Nonparametric estimation: parametric view.

[10]
Stefan Steinerberger.
Wasserstein distance, fourier series and applications.
arXiv:1803.08011, 2018.

[11]
Anja Sturm Thomas Rippl, Axel Munk.
Limit laws of the empirical wasserstein distance: Gaussian
arXiv:1507.04090v2, 2015.