12

# Information Geometry of Wasserstein Statistics
on Shapes and Affine Deformations

Shun-ichi Amari
Teikyo University, Advanced Comprehensive Research Organization,
RIKEN Center for Brain Science  Takeru Matsuda
The University of Tokyo,
RIKEN Center for Brain Science

###### Abstract

Information geometry and Wasserstein geometry are two main structures
introduced in a manifold of probability distributions, and they capture its
different characteristics. We study characteristics of Wasserstein geometry in
the framework of Li and Zhao (2023) for the affine deformation statistical
model, which is a multi-dimensional generalization of the location-scale
model. We compare merits and demerits of estimators based on information
geometry and Wasserstein geometry. The shape of a probability distribution and
its affine deformation are separated in the Wasserstein geometry, showing its
robustness against the waveform perturbation in exchange for the loss in
Fisher efficiency. We show that the Wasserstein estimator is the moment
estimator in the case of the elliptically symmetric affine deformation model.
It coincides with the information-geometrical estimator (maximum-likelihood
estimator) when and only when the waveform is Gaussian. The role of the
Wasserstein efficiency is elucidated in terms of robustness against waveform
change.

##  1 Introduction

We study statistics based on probability distribution patterns $p({\bm{x}})$
over ${\bm{x}}\in X={\bm{R}}^{d}$, by using both information geometry (see
Amari, 2016; Ay et al., 2017, etc) and Wasserstein geometry (see Villani,
2003; Peyré and Cuturi, 2019; Santambrogio, 2015, among many others). Here,
$p({\bm{x}})$ is a probability distribution on $X={\bm{R}}^{d}$. When $d=2$,
$p({\bm{x}})$ is regarded as a visual pattern on ${\bm{R}}^{2}$.

There are lots of applications of Wasserstein geometry to statistics (see,
e.g., Amari and Matsuda, 2022; Bernton et al., 2019; Yatracos, 2022; Bassetti
et al., 2006; Matsuda and Strawderman, 2021; Imaizumi et al., 2022; Li and
Montufar, 2020; Chen et al., 2021, and others), machine learning (see, e.g.,
Arjovsky et al., 2017; Fronger et al., 2015; Wang and Li, 2020; Peyré and
Cuturi, 2019; Montavon et al., 2015, among many others) and statistical
physics (Ito, 2023). We have a good review paper and a book, Panaretos and
Zemel (2019, 2022), which include lots of references. However, applications to
statistical inference look still premature. For example, the statistical
efficiency of the Wasserstein estimator is studied only in the one-dimensional
location-scale model (Amari and Matsuda, 2022). We give characterization of
the Wasserstein estimator from the point of view of the robustness for changes
of the shape of probability distributions. We further focus on the efficiency
of the Wasserstein estimator for the affine deformation statistical model
(deformation model in short), where deformation parameters and the waveform of
a probability distribution are separated.

The affine deformation statistical model $p({\bm{x}},{\bm{\theta}})$ is
generated from a standard shape distribution $f({\bm{z}})$ satisfying

$\displaystyle\int f({\bm{z}})d{\bm{z}}$ $\displaystyle=$ $\displaystyle 1,$
(1) $\displaystyle\int{\bm{z}}f({\bm{z}})d{\bm{z}}$ $\displaystyle=$
$\displaystyle 0,$ (2)
$\displaystyle\int{\bm{z}}{\bm{z}}^{\top}f({\bm{z}})d{\bm{z}}$
$\displaystyle=$ $\displaystyle I,$ (3)

where $I$ is the identity matrix. The deformation parameter consists of
${\bm{\theta}}=\left({\bm{\mu}},\Lambda\right)\in\Theta$ such that
${\bm{\mu}}$ is a vector specifying translation of the location and $\Lambda$
is a non-singular matrix representing scale changes and rotations of
${\bm{x}}$. Then, $p({\bm{x}},{\bm{\theta}})$ is written as

$p({\bm{x}},{\bm{\theta}})=|\Lambda|f\left(\Lambda({\bm{x}}-{\bm{\mu}})\right).$
(4)

Given a standard $f$, we have a statistical model parameterized by
${\bm{\theta}}$,

$M_{f}=\left\\{p({\bm{x}},{\bm{\theta}})\right\\}.$ (5)

Geometrically, it forms a finite-dimensional statistical manifold, where
${\bm{\theta}}$ plays the role of a coordinate system. The deformation model
is a generalization of the location-scale model. Note that this model is often
called the location-scatter model in several fields such as statistics and
signal processing (Tyler, 1987; Ollila and Tyler, 2014).

Let $T_{\bm{\theta}}$ denote the affine deformation from ${\bm{x}}$ to
${\bm{z}}$ given by

${\bm{z}}=T_{\bm{\theta}}{\bm{x}}=\Lambda({\bm{x}}-{\bm{\mu}}).$ (6)

This may be regarded as deformation of shape $f$ to
$\tilde{T}_{\bm{\theta}}f$,

$(\tilde{T}_{\bm{\theta}}f)({\bm{x}})=f\left(T_{\bm{\theta}}{\bm{x}}\right),$
(7)

where $\tilde{T}_{\bm{\theta}}$ is an operator to change a standard waveform
$f$ to another waveform $\tilde{f}=\tilde{T}_{\bm{\theta}}f$ defined by (7).

Let ${\cal{F}}=\\{q({\bm{x}})\\}$ be the space of all smooth positive
probability density functions that have mean and covariance. Let
${\cal{F}}_{S}=\\{f({\bm{z}})\\}$ be its subspace consisting of all the
standard distributions $f({\bm{z}})$ satisfying (1), (2) and (3). Then, any
$q({\bm{x}})\in{\cal{F}}$ is uniquely written in the form

$q({\bm{x}})=|\Lambda|f\left(\Lambda({\bm{x}}-{\bm{\mu}})\right)$ (8)

for $f\in{\cal{F}}_{S}$ and ${\bm{\theta}}=({\bm{\mu}},\Lambda)\in\Theta$.
Hence, ${\cal{F}}$ is the direct product of ${\cal{F}}_{S}$ and $\Theta$:
${\cal{F}}={\cal{F}}_{S}\times\Theta$. See Figure 1.

Figure 1: Decomposition of $\cal{F}$

It is interesting to see that ${\cal{F}}_{S}$ is not a manifold. It includes
lots of singularities. Let $O$ be an orthogonal transformation. Then,
$f({\bm{x}})$ and $(Of)({\bm{x}})=f(O{\bm{x}})$ give different shape functions
in general. However, when $f$ is rotationally invariant, the shapes of $f$ and
$Of$ are identical, $f=Of$. The isotropic Gaussian shape is such an example,
and all $Of$ reduces to one point in such a case. It is also the case when $f$
is invariant under a subgroup of the orthogonal group. See Figure 2 for the
shape of ${\cal{F}}_{S}$. It is a covex set, because, for
$f,g\in{\cal{F}}_{S}$, $\alpha f+(1-\alpha)g$ also belongs to ${\cal{F}}_{S}$
for $0\leq\alpha\leq 1$.

Figure 2: Singular structure of ${\cal{F}}_{S}$

Geometry of a manifold of probability distributions has so far been studied by
information geometry and Wasserstein geometry. The two geometries capture
different aspects of a manifold of probability distributions. We use a
divergence measure to explain this. Let $D_{F}[p({\bm{x}}),q({\bm{x}})]$ and
$D_{W}[p({\bm{x}}),q({\bm{x}})]$ be two divergence measures between
distributions $p({\bm{x}})$ and $q({\bm{x}})$, where subscripts $F$ and $W$
represent Fisher-based information geometry and Wasserstein geometry,
respectively. Information geometry uses an invariant divergence $D_{F}$,
typically the Kullback–Leibler divergence. Wasserstein divergence $D_{W}$ is
defined by the cost of transporting masses distributed in form $p({\bm{x}})$
to another $q({\bm{x}})$. Roughly speaking, $D_{F}$ measures the vertical
differences of $p({\bm{x}})$ and $q({\bm{x}})$, for example, represented by
their log-ratio $\log(p({\bm{x}})/q({\bm{x}}))$, at ${\bm{x}}$, whereas
$D_{W}$ measures the horizontal differences of $p({\bm{x}})$ and $q({\bm{x}})$
which corresponds to the transportation cost from $p({\bm{x}})$ to
$q({\bm{x}})$. See Figure 3.

Figure 3: (a) $F$-divergence. (b) $W$-divergence.

Information geometry is constructed based on the invariance principle of
Chentsov (Chentsov, 1982) such that $D_{F}[p({\bm{x}}),q({\bm{x}})]$ is
invariant under invertible transformations of the coordinates ${\bm{x}}$ of
the sample space $X$. This implies that the divergence does not depend on the
coordinate system of $X$. We then have a unique Riemannian metric, which is
Fisher–Rao metric, and also a dual pair of affine connections (Amari and
Nagaoka, 2007). This is useful not only for analyzing the performances of
statistical inference but also for vision analysis, machine learning,
statistical physics, and many others (see Amari, 2016).

Wasserstein geometry has an old origin, proposed by G. Monge in 1781 as a
problem of transporting mass distributed in the form $p({\bm{x}})$ to another
$q({\bm{x}})$ such that the total transportation cost is minimized. It depends
on the transportation cost $c({\bm{x}},{\bm{y}})$ between two locations
${\bm{x}},{\bm{y}}\in X$. The cost is usually a function of the Euclidean
distance between ${\bm{x}}$ and ${\bm{y}}$. We use the square of the distance
as a cost function, which gives $L^{2}$-Wasserstein geometry. This Wasserstein
geometry directly depends on the Euclidean distance of $X={\bm{R}}^{d}$.
Therefore, it is responsible for the metric structure of $X$ and is useful for
problems that intrinsically depend on the structure of $X$, such as the
transportation problem, non-equilibrium statistical physics, pattern analysis,
machine learning and many others.

It is natural to search for the relation between the two geometries. There are
a number of such trials, including Amari et al. (2018, 2019); Khan et al.
(2022); Rankin and Wong (2023); Ito (2023) and others. Among them, Li and Zhao
(2023) gave a unified framework for the two geometries. The present article is
based on their framework and focuses on the affine deformation model, for
which the standard waveform $f$ and the deformation parameter ${\bm{\theta}}$
are separated.

Li and Zhao (2023) introduced the $W$-score function in parallel to the Fisher
score function, defining two estimators $\hat{\bm{\theta}}_{F}$ and
$\hat{\bm{\theta}}_{W}$ thereby. The former is the maximum likelihood
estimator that maximizes the log likelihood. This is the one that minimizes an
invariant divergence from the empirical distribution $\hat{p}({\bm{x}})$ to
parametrized model $M_{f}$, where the empirical distribution is given based on
$n$ independent and identical observations ${\bm{x}}_{1},\cdots,{\bm{x}}_{n}$
as

$\hat{p}({\bm{x}})=\frac{1}{n}\sum_{i}\delta\left({\bm{x}}-{\bm{x}}_{i}\right),$
(9)

and $\delta({\bm{x}})$ is the delta function. The latter $W$-estimator
$\hat{\bm{\theta}}_{W}$ is asymptotically equivalent to the minimizer of the
$W$-divergence between the empirical distribution and model $M_{f}$ (see
Section 4). Li and Zhao (2023) further defined the $F$-efficiency and
$W$-efficiency of an estimator $\hat{\bm{\theta}}$ given a statistical model
$M=\left\\{p({\bm{x}},{\bm{\theta}})\right\\}$, proving the Cramér-Rao type
inequalities. We apply their theory to analyze the effects of the shape $f$ on
the efficiencies of the two types of estimators.

The present paper is organized as follows. In section 2, we introduce two
divergences between distributions, one based on the invariance principle and
the other based on the transportation cost. The divergences give two
Riemannian structures in the space ${\cal{F}}$ of probability distributions
$p({\bm{x}})$ over $X={\bm{R}}^{d}$. A regular statistical model
$M=\left\\{p({\bm{x}},{\bm{\theta}})\right\\}$ parameterized by
${\bm{\theta}}$ is a finite-dimensional submanifold embedded in ${\cal{F}}$.
In section 3, we define the $F$\- and $W$-score functions following Li and
Zhao (2023). The Riemannian structure of the tangent space of probability
distributions is pulled-back to the model submanifold, giving both the
Riemannian metrics and score functions to $M$. We define the $F$\- and
$W$-estimators $\hat{\bm{\theta}}_{F}$ and $\hat{\bm{\theta}}_{W}$ by using
the $F$\- and $W$-score functions, respectively. Section 4 defines the affine
deformation statistical model. Section 5 studies the elliptically symmetric
affine deformation model $M_{f}$, where $f$ is a spehrically symmetric
standard form. For this model, we show that the $W$-score functions are
quadratic functions of ${\bm{x}}$. Hence, it is proved that
$\hat{\bm{\theta}}_{W}$ is a moment estimator. We also show that $M_{f}$ and
${\cal{F}}_{S}$ are orthogonal in the $W$-geometry, implying the separation of
the waveform and deformation. In Section 6, we elucidate the role of
$W$-efficiency from the point of view of robustness to a change in the
waveform $f$ due to observation noise. In Section 7, we prove that the
Gaussian shape is a unique model in which $F$-estimator and $W$-estimator
coincide, with $\hat{\bm{\theta}}_{W}$ satisfying the $F$-efficiency and
$\hat{\bm{\theta}}_{F}$ satisfying $W$-efficiency. Section 8 briefly
summarizes the paper and mentions future work.

##  2 Riemannian structures in the space of probability densities on
${\bm{R}}^{d}$

We consider the space ${\cal{F}}=\left\\{p({\bm{x}})\right\\}$ of all smooth
positive probability density functions on ${\bm{R}}^{d}$, having the mean and
covariance. Later, we may relax the conditions of positivity and smoothness,
when we discuss a parametric model, in particular the deformation model. We
define a divergence function $D[p({\bm{x}}),q({\bm{x}})]$, which represents
the degree of difference between $p({\bm{x}})$ and $q({\bm{x}})$. The square
of the distance between $p({\bm{x}})$ and $q({\bm{x}})$ plays this role, but a
divergence does not necessarily need to be symmetric with respect to
$p({\bm{x}})$ and $q({\bm{x}})$. A divergence function satisfies the following
conditions:

  1. 1.

$D[p({\bm{x}}),q({\bm{x}})]\geq 0$ and the equality holds if and only if
$p({\bm{x}})=q({\bm{x}})$.

  2. 2.

Let $\delta p({\bm{x}})$ be an infinitesimally small deviation of
$p({\bm{x}})$. Then, $D[p({\bm{x}}),p({\bm{x}})+\delta p({\bm{x}})]$ is
approximated by a positive quadratic functional of $\delta p({\bm{x}})$.

A divergence is said to be invariant if

$D\left[p({\bm{x}}),q({\bm{x}})\right]=D\left[\tilde{p}({\bm{y}}),\tilde{q}({\bm{y}})\right]$
(10)

holds for every smooth reversible transformation ${\bm{k}}$ of the coordinates
from ${\bm{x}}\in{\bm{R}}^{d}$ to ${\bm{y}}={\bm{k}}({\bm{x}})$, where

$\tilde{p}({\bm{y}})=\left|\frac{\partial{\bm{x}}}{\partial{\bm{y}}}\right|p({\bm{x}}).$
(11)

A typical invariant divergence is the $\alpha$-divergence $(\alpha\neq\pm 1)$
defined by

$\displaystyle D_{\alpha}[p,q]$ $\displaystyle=$
$\displaystyle\frac{4}{1-\alpha^{2}}\left(1-\int
p({\bm{x}})^{(1+\alpha)/2}q({\bm{x}})^{(1-\alpha)/2}d{\bm{x}}\right)$ (12)

for $\alpha\neq\pm 1$. For $\alpha=1$, we define $D_{1}[p,q]$ by the
Kullback–Leibler divergence

$D_{1}[p,q]=\int p({\bm{x}})\log\frac{p({\bm{x}})}{q({\bm{x}})}d{\bm{x}}.$
(13)

For $\alpha=-1$, we define $D_{-1}[p:q]=D_{1}[q:p]$. The case $\alpha=0$ is
equivalent to the Hellinger divergence

$H^{2}[p,q]=\frac{1}{2}\int\left(\sqrt{p({\bm{x}})}-\sqrt{q({\bm{x}})}\right)^{2}d{\bm{x}}.$
(14)

A characterization of the $\alpha$-divergence is given in Amari (2016). The
$\alpha$-divergence gives information-geometric structure to $\cal{F}$.

Another divergence is the Wasserstein divergence. Let us transport masses
piled in the form $p({\bm{x}})$ to another $q({\bm{x}})$. To this end, we need
to move some mass at ${\bm{x}}$ to another position ${\bm{y}}$. Let
$\pi({\bm{x}},{\bm{y}})$ be a stochastic matrix, showing the probability of
mass at ${\bm{x}}$ to be transported to ${\bm{y}}$. We call $\pi$ a
transportation plan when it satisfies the following terminal conditions

$\displaystyle\int\pi({\bm{x}},{\bm{y}})d{\bm{y}}$ $\displaystyle=$
$\displaystyle p({\bm{x}}),$ (15)
$\displaystyle\int\pi({\bm{x}},{\bm{y}})d{\bm{x}}$ $\displaystyle=$
$\displaystyle q({\bm{y}}).$ (16)

Let $c({\bm{x}},{\bm{y}})$ be the cost of transporting a unit of mass from
${\bm{x}}$ to ${\bm{y}}$. Then, the Wasserstein divergence
$D_{W}[p({\bm{x}}),q({\bm{x}})]$ is the minimum transporting cost from
$p({\bm{x}})$ to $q({\bm{x}})$. By using stochastic plan
$\pi({\bm{x}},{\bm{y}})$, the Wasserstein divergence between $p({\bm{x}})$ and
$q({\bm{x}})$ is given by

$D_{W}\left[p({\bm{x}}),q({\bm{x}})\right]={\mathop{\inf}_{\pi}}\int
c({\bm{x}},{\bm{y}})\pi({\bm{x}},{\bm{y}})d{\bm{x}}d{\bm{y}},$ (17)

where infimum is taken over all stochastic plans $\pi$ satisfying (15) and
(16). When the cost is the square of the Euclidean distance

$c({\bm{x}},{\bm{y}})=\|{\bm{x}}-{\bm{y}}\|^{2},$ (18)

we call $D_{W}$ the $L^{2}$-Wasserstein divergence. We focus on this
divergence in the following. Note that the $L^{2}$-Wasserstein divergence is
the square of the $L^{2}$-Wasserstein distance.

The dynamic formulation of the optimal transport problem proposed by Brenier
(1999) and developed further by Benamou and Brenier (2000) is useful. Let
$\rho({\bm{x}},t)$ be a family of probability distributions parameterized by
$t$. It represents the time course $\rho({\bm{x}},t)$ of transporting
$p({\bm{x}})$ to $q({\bm{x}})$, satisfying

$\rho({\bm{x}},0)=p({\bm{x}}),\quad\rho({\bm{x}},1)=q({\bm{x}}).$ (19)

We introduce potential $\Phi({\bm{x}},t)$ such that its gradient
$\nabla_{\bm{x}}\Phi({\bm{x}},t)$ represents the velocity

${\bm{v}}({\bm{x}},t)=\nabla_{\bm{x}}\Phi({\bm{x}},t)$ (20)

of mass flow at ${\bm{x}}$ and $t$ in the dynamic plan. Then,
$\Phi({\bm{x}},t)$ satisfies the following equation of continuity

$\partial_{t}\rho({\bm{x}},t)+\nabla_{\bm{x}}\cdot\left\\{\rho({\bm{x}},t)\nabla_{\bm{x}}\Phi({\bm{x}},t)\right\\}=0.$
(21)

The Wasserstein divergence is written in the dynamic formulation as

$D_{W}\left[p({\bm{x}}),q({\bm{x}})\right]={\mathop{\inf}_{\Phi}}\int^{1}_{0}\int\|\nabla_{\bm{x}}\Phi({\bm{x}},t)\|^{2}\rho({\bm{x}},t)d{\bm{x}}dt.$
(22)

We introduce a Riemannian structure to $\cal{F}$ by the Taylor expansion of
$D[p,p+\delta p]$. The Riemannian metric $g$ is an operator to give the
squared magnitude $ds^{2}$ of an infinitensimal deviation $\delta p({\bm{x}})$
in the tangent space of $\cal{F}$, for example, by

$ds^{2}=\int\delta p({\bm{x}})g({\bm{x}},{\bm{y}})\delta
p({\bm{y}})d{\bm{x}}d{\bm{y}}.$ (23)

In the case of the invariant divergence, we have

$g_{F}({\bm{x}},{\bm{y}})=\frac{\delta{({\bm{x}}-{\bm{y}})}}{p({\bm{x}})},$
(24)

where $\delta({\bm{x}})$ is the delta function and $g_{F}$ is a positive
integral operator. In the case of the $L^{2}$-Wasserstein divergence, for
infinitesimally small $\delta t$, the change of distribution from
$\rho({\bm{x}},0)=p({\bm{x}})$ at $t=0$ to $\rho({\bm{x}},\delta
t)=p({\bm{x}})+\delta p({\bm{x}})$ at $t=\delta t$ satisfies

$\delta p({\bm{x}})=-\Delta_{p}\Phi({\bm{x}})dt,$ (25)

where $\Delta_{p}$ is the $p$-Laplacian defined by

$\Delta_{p}\Phi=\nabla_{\bm{x}}\cdot\left(p({\bm{x}})\nabla_{\bm{x}}\Phi({\bm{x}})\right).$
(26)

The $L^{2}$-Wasserstein divergence is given in the quadratic form of $\delta
p$ as

$D_{W}[p,p+\delta p]=-\int\delta p({\bm{x}})\Delta^{-1}_{\rho}\delta
p({\bm{x}})d{\bm{x}},$ (27)

giving the Riemannian metric $g_{W}$ which is an operator represented by

$g_{W}(p)=-\Delta^{-1}_{p}.$ (28)

This is Otto’s Riemannian metric (Otto, 2001). See Li and Zhao (2023) for
details.

##  3 Score functions and estimators in parametric model

We consider a regular statistical model
$M=\left\\{p({\bm{x}},{\bm{\theta}})\right\\}$ parameterized by
$m$-dimensional vector ${\bm{\theta}}$. The tangent space of $M$ at
${\bm{\theta}}$ is spanned by $\partial_{i}p({\bm{x}},{\bm{\theta}})$ for
$i=1,\cdots,m$, such that

$\delta p({\bm{x}})=\partial_{i}p({\bm{x}},{\bm{\theta}})d\theta^{i},$ (29)

where

$\partial_{i}p({\bm{x}},{\bm{\theta}})=\frac{\partial}{\partial\theta^{i}}p({\bm{x}},{\bm{\theta}}).$
(30)

Hereafter, the summation convention is used, that is, all indices appearing
twice, once as upper and the other as lower indices, e.g. $i$’s in (29), are
summed up.

Let us define $S_{i}({\bm{x}},{\bm{\theta}})$ from the basis functions
$\partial_{i}p({\bm{x}},{\bm{\theta}})$ of the tangent space of $M$ for
$i=1,\cdots,m$ by applying the Riemannian metric operator $g$:

$S_{i}({\bm{x}},{\bm{\theta}})=g\circ\partial_{i}p({\bm{x}},{\bm{\theta}}).$
(31)

We call them score functions following the tradition of statistics. In the
case of invariant Fisher geometry, the score functions are

$\displaystyle S^{F}_{i}({\bm{x}},{\bm{\theta}})$ $\displaystyle=$
$\displaystyle\frac{\partial_{i}p}{p({\bm{x}},{\bm{\theta}})}=\partial_{i}l({\bm{x}},{\bm{\theta}}),$
(32) $\displaystyle l({\bm{x}},{\bm{\theta}})$ $\displaystyle=$
$\displaystyle\log p({\bm{x}},{\bm{\theta}}),$ (33)

which is the derivative of log-likelihood. In Wasserstein geometry, we have

$S^{W}_{i}({\bm{x}},{\bm{\theta}})=-\Delta_{p}^{-1}\partial_{i}p({\bm{x}},{\bm{\theta}}).$
(34)

By using the identity

$\int
a({\bm{x}})\Delta_{p}b({\bm{x}})d{\bm{x}}=-\int\left(\nabla_{\bm{x}}a\cdot\nabla_{\bm{x}}b\right)p({\bm{x}})d{\bm{x}},$
(35)

we see that $S^{W}_{i}({\bm{x}},{\bm{\theta}})$ are solutions of the following
Poisson equations,

$\nabla_{\bm{x}}\log
p({\bm{x}},{\bm{\theta}})\cdot\nabla_{\bm{x}}S^{W}_{i}({\bm{x}},{\bm{\theta}})+\Delta_{\bm{x}}S^{W}_{i}({\bm{x}},{\bm{\theta}})+\frac{\partial}{\partial\theta_{i}}\log
p({\bm{x}},{\bm{\theta}})=0.$ (36)

For infinitesimal $\delta$, the map
$\bm{x}\mapsto\bm{x}+\delta\nabla_{\bm{x}}S_{i}^{{W}}(\bm{x},\bm{\theta})$ is
the optimal transport map from $p(\bm{x},\bm{\theta})$ to
$p(\bm{x},\bm{\theta}+\delta\bm{e}_{i})$ with transportation cost

$D_{W}(p(\bm{x},\bm{\theta}),p(\bm{x},\bm{\theta}+\delta\bm{e}_{i}))=\int\|\delta\nabla_{\bm{x}}S_{i}^{{W}}(\bm{x},\bm{\theta})\|^{2}p(\bm{x},\bm{\theta}){\rm
d}x,$ (37)

where $\bm{e}_{i}$ is the $i$-th standard unit vector. In order to eliminate
the indefiniteness due to the integral constant, we pose an additional
condition

${\rm{E}}\left[S_{i}({\bm{x}},{\bm{\theta}})\right]=0.$ (38)

This is automatically satisfied in the Fisherian case.

The Riemannian metric tensor $g_{ij}({\bm{\theta}})$ is pulled-back from $g$
in ${\cal{F}}$ to $M$, and is derived in terms of the score functions as

$g_{ij}({\bm{\theta}})=\langle S_{i},g^{-1}S_{j}\rangle,$ (39)

where

$\langle a({\bm{x}}),b({\bm{x}})\rangle=\int a({\bm{x}})b({\bm{x}})d{\bm{x}}.$
(40)

In the Fisherian case,

$g^{F}_{ij}({\bm{\theta}})={\rm{E}}\left[\partial_{i}l\partial_{j}l\right]=\int
p(x,{\bm{\theta}})\partial_{i}l({\bm{x}},{\bm{\theta}})\partial_{j}l({\bm{x}},{\bm{\theta}})d{\bm{x}}.$
(41)

In the Wasserstein case,

$g^{W}_{ij}({\bm{\theta}})=\int
p({\bm{x}},{\bm{\theta}})\nabla_{\bm{x}}S^{W}_{i}({\bm{x}},{\bm{\theta}})^{\top}\nabla_{\bm{x}}S^{W}_{j}({\bm{x}},{\bm{\theta}})d{\bm{x}}={\rm
E}[\nabla_{\bm{x}}S^{W}_{i}({\bm{x}},{\bm{\theta}})^{\top}\nabla_{\bm{x}}S^{W}_{j}({\bm{x}},{\bm{\theta}})],$
(42)

where identity (35) is used.

The score functions $S_{i}({\bm{x}},{\bm{\theta}})$ give a set of estimating
functions from (38), which are used to obtain an estimator
$\hat{\bm{\theta}}$. Let $\hat{p}_{\rm{emp}}({\bm{x}})$ be the empirical
distribution given by

$\hat{p}_{\rm{emp}}({\bm{x}})=\frac{1}{n}\sum^{n}_{j=1}\delta\left({\bm{x}}-{\bm{x}}_{j}\right),$
(43)

where ${\bm{x}}_{1},\cdots,{\bm{x}}_{n}$ are $n$ independent observations.
Then, replacing expectation ${\rm{E}}$ in (38) by the expectation with respect
to the empirical distribution, we have estimating equations,

${\rm{E}}_{\rm{emp}}\left[S_{i}({\bm{x}},{\bm{\theta}})\right]=\frac{1}{n}\sum^{n}_{j=1}S_{i}\left({\bm{x}}_{j},{\bm{\theta}}\right)=0,\quad
i=1,\cdots,m.$ (44)

It is known that the solution $\hat{\bm{\theta}}$ gives a consistent estimator
for large $n$. Roughly speaking, $\hat{\bm{\theta}}$ is the projection of
$\hat{p}_{\rm{emp}}$ to the model $M$ with respect to the metric $g$ (see
Figure 4). It is the solution of

$\langle\hat{p}_{\rm{emp}}({\bm{x}}),S_{i}({\bm{x}},{\bm{\theta}})\rangle=0,$
(45)

giving a consistent estimator $\hat{\bm{\theta}}$. A consistent estimator is
Fisher efficient when the projection is orthogonal with respect to Fisher-Rao
metric (Amari and Nagaoka, 2007).

Figure 4: Projection of $\hat{p}_{\rm{emp}}$ to $M$.

Estimator $\hat{\bm{\theta}}_{F}$ (the invariant Fisherian case) is the
maximum likelihood estimator, that maximizes the likelihood. Cramér-Rao
theorem gives a matrix inequality for any unbiased estimator
$\hat{\bm{\theta}}$,

${\rm{Cov}}\left[\hat{\bm{\theta}}-{\bm{\theta}}\right]\succeq\frac{1}{n}g^{-1}_{F}({\bm{\theta}}),$
(46)

where ${\rm{Cov}}[\cdot]$ is the covariance matrix and $\succeq$ denotes the
matrix order defined by the positive definiteness. The maximum likelihood
estimator $\hat{\bm{\theta}}_{F}$ satisfies

${\rm{Cov}}\left[\hat{\bm{\theta}}_{F}-{\bm{\theta}}\right]\approx\frac{1}{n}g^{-1}_{F}({\bm{\theta}})$
(47)

asymptotically. Hence, it minimizes the error covariance matrix and the
minimized error covariance is given asymptotically by the inverse of the
Fisher metric tensor $g_{F}$ divided by $n$. Such a property is called the
Fisher efficiency.

In parallel to the Fisherian case, we need to study the characteristics of the
Wasserstein estimator $\hat{\bm{\theta}}_{W}$ in the following. In the case of
the one-dimensional location-scale model, the Wasserstein estimator is
asymptotically equivalent to the estimator obtained by minimizing the
Wasserstein divergence (transportation cost) from the empirical distribution
$\hat{p}_{\rm{emp}}$ to model $M$,

$\hat{\bm{\theta}}_{W}={\mathop{\arg\min}_{{\bm{\theta}}}}\;D_{W}\left[\hat{p}_{\rm{emp}},p({\bm{x}},{\bm{\theta}})\right].$
(48)

See the end of Section 5. The properties of $\hat{\bm{\theta}}_{W}$ were
studied in detail by Amari and Matsuda (2022) in the case of the one-
dimensional location-scale model.

##  4 Affine deformation model

Now, we focus on the affine deformation model. Let $f({\bm{z}})$ be a standard
probability density function satisfying (1), (2), and (3). To define $M_{f}$,
we use affine deformation of ${\bm{x}}$ to ${\bm{z}}$ by

${\bm{z}}=\Lambda({\bm{x}}-{\bm{\mu}}),$ (49)

where ${\bm{\mu}}$ is a vector representing shift of location and $\Lambda$ is
a non-singular matrix. Hence, ${\bm{\theta}}=({\bm{\mu}},\Lambda)$ is
$m=d^{2}+d$ dimensional. The model $M_{f}$ is defined from

$p({\bm{x}},{\bm{\theta}})d{\bm{x}}=f({\bm{z}})d{\bm{z}},$ (50)

that is

$p({\bm{x}},{\bm{\theta}})=|\Lambda|f(\Lambda({\bm{x}}-{\bm{\mu}})),$ (51)

satisfying

$\displaystyle\int{\bm{x}}p({\bm{x}},{\bm{\theta}})d{\bm{x}}$ $\displaystyle=$
$\displaystyle{\bm{\mu}},$ (52)
$\displaystyle\int{\bm{x}}{\bm{x}}^{\top}p({\bm{x}},{\bm{\theta}})d{\bm{x}}$
$\displaystyle=$ $\displaystyle\Lambda^{-2}+{\bm{\mu}}{\bm{\mu}}^{\top}.$ (53)

This is a generalization of the location-scale model, which is simply obtained
by putting $\Lambda=(1/\sigma)I$, with $\sigma$ being the scale factor. It
should be noted that $\Lambda$ is decomposed as $\Lambda=UDO$, where $U$ and
$O$ are orthogonal matrices and $D$ is a positive diagonal matrix. In the
following, we denote the log probability of standard shape $f$ by

$l({\bm{z}})=\log f({\bm{z}}).$ (54)

As we discussed in Introduction, the set of all standard shape functions
${\cal{F}}_{S}=\\{f\\}$ does not form a manifold but has an interesting
topological structure due to the rotational invariance for some $f$. For each
standard shape function $f\in{\cal{F}}_{S}$, an affine deformation model
$M_{f}$ parameterized by ${\bm{\theta}}=({\bm{\mu}},\Lambda)$ is attached.
Thus, ${\cal{F}}$ is decomposed into the direct product of ${\cal{F}}_{S}$ and
$M_{f}$,

${\cal{F}}={\cal{F}}_{S}\times M_{f}.$ (55)

For any $f$, $M_{f}$ has cone structure parameterized by $({\bm{\mu}},D,U,O)$,
where $\Lambda=UDO$ and $D$ is a diagonal matrix with diagonal elements
$d_{i}>0$. Thus, $D$ can be identified with a vector in the open positive
quadrant $\bm{R}_{+}^{d}$ of ${\bm{R}}^{d}$, which has the cone structure.
Since ${\bm{\mu}}\in{\bm{R}}^{d}$, and $U,O\in{\cal{O}}(d)$, we have the
decomposition

$M_{f}={\bm{R}}^{d}\times\bm{R}_{+}^{d}\times{\cal{O}}(d)\times{\cal{O}}(d).$
(56)

See Takatsu and Yokota (2012) for the cone structure of ${\cal{F}}$. When $f$
is Gaussian, its structure is studied in detail by Takatsu (2011).

When $p({\bm{x}})$ belongs to $M_{f}$, the waveform of $p({\bm{x}})$ is said
to be equivalent to that of $f$. $M_{f}$ consists of the distributions of all
equivalent waveforms. All ellipsoidal shapes are equivalent to a spherical
shape. A family of special parallel-piped shapes are equivalent to a cubic
form (see Fig. 5). Therefore, our model is useful for separating the effect of
the shape from location and affine deformation.

Figure 5: Equivalent shapes.

We may consider subclasses of the transformation model. One simple example is
the location model, in which $\Lambda$ is fixed to the identity matrix $I$. A
stronger theorem is known in such a simple model (Givens and Shortt, 1984). In
our context, it can be expressed as follows.

######  Proposition 1.

Wasserstein geometry gives an orthogonal decomposition of the shape and
locations,

$D^{W}\left[f_{1}\left({\bm{x}}-{\bm{\mu}}_{1}\right),f_{2}\left({\bm{x}}-{\bm{\mu}}_{2}\right)\right]=D^{W}\left[f_{1}({\bm{x}}),f_{2}({\bm{x}})\right]+\|{\bm{\mu}}_{1}-{\bm{\mu}}_{2}\|^{2}.$
(57)

##  5 Elliptically symmetric deformation model

Here, we focus on deformation models that are elliptically symmetric:

$p({\bm{x}},{\bm{\theta}})=|\Lambda|g(\|\Lambda({\bm{x}}-{\bm{\mu}})\|),$ (58)

where $f(\bm{z})=g(\|\bm{z}\|)$ satisfies the standard density conditions (1),
(2), and (3). Note that $f(\bm{z})=f(UDO\bm{x})=g(\|DO\bm{x}\|)$ does not
depend on $U$ and thus the parameter is reduced to $(\bm{\mu},D,O)$ in this
case.

First, we consider the $F$-estimator $\hat{\bm{\theta}}_{F}$ (maximum
likelihood estimator). The log-likelihood is given by

$\log p({\bm{x}},{\bm{\theta}})=\log|\Lambda|+\log
g(\|\Lambda({\bm{x}}-{\bm{\mu}})\|).$ (59)

When there are $n$ observations ${\bm{x}}_{1},\cdots,{\bm{x}}_{n}$, summation
is taken over them so that we have the likelihood equations

$\sum^{n}_{j=1}\partial_{\bm{\theta}}\log
p\left({\bm{x}}_{j},{\bm{\theta}}\right)=0.$ (60)

The solution $\hat{\bm{\theta}}_{F}$ strongly depends of the shape $g$.

Contrary to this, the $W$-estimator $\hat{\bm{\theta}}_{W}$ does not depend on
the shape $g$ as follows.

######  Lemma 1.

$\displaystyle\nabla_{x}\left(\frac{1}{2}x^{\top}Ax+b^{\top}x+c\right)=\frac{A+A^{\top}}{2}x+b,\quad\Delta_{x}\left(\frac{1}{2}x^{\top}Ax+b^{\top}x+c\right)={\rm
tr}(A).$

###### Proof.

Straightforward calculation. Note that $(A+A^{\top})/2\neq A$ when $A$ is not
symmetric. ∎

######  Lemma 2.

Let $A,B\in\mathbb{R}^{d\times d}$ be symmetric matrices. If $A$ is positive
semidefinite, then the Sylvester equation $AX+XA=B$ has the unique solution
$X$, which satisfies $X^{\top}=X$ and ${\rm tr}(X)={\rm tr}(A^{-1}B)/2$.

###### Proof.

From the positive semidefiniteness of $A$, the spectra of $A$ and $-A$ are
disjoint. Thus, from Theorem VII.2.1 of Bhatia (1997), the Sylvester equation
$AX+XA=B$ has a unique solution.

Let $X$ be the solution of the Sylvester equation. From $A^{\top}=A$, we have
$AX^{\top}+X^{\top}A=(AX+XA)^{\top}=B^{\top}=B$, which means that $X^{\top}$
is also a solution of the Sylvester equation. Since the solution is unique, it
implies $X^{\top}=X$. Also, from the positive semidefiniteness of $A$ and
$AX+XA=B$, we have $X+A^{-1}XA=A^{-1}B$. Taking the trace and using ${\rm
tr}(A^{-1}XA)={\rm tr}(X)$, we obtain ${\rm tr}X={\rm tr}(A^{-1}B)/2$. ∎

######  Theorem 1.

For the elliptically symmetric deformation model (58), the Wasserstein score
functions are quadratic. Specifically, the Wasserstein score function for
$\mu_{i}$ is

$\displaystyle S_{\mu_{i}}^{{W}}(\bm{x},\bm{\theta})=x_{i}-\mu_{i},$

and the Wasserstein score function for $\Lambda_{ij}$ is

$\displaystyle
S_{\Lambda_{ij}}^{W}(\bm{x},\bm{\theta})=\frac{1}{2}\bm{x}^{\top}A\bm{x}+b^{\top}\bm{x}-{\rm
E}_{\bm{\theta}}\left[\frac{1}{2}\bm{x}^{\top}A\bm{x}+b^{\top}\bm{x}\right],$

where $A$ is the unique solution of the Sylvester equation
$\Lambda^{2}A+A\Lambda^{2}=-L$ and $b=-A\mu$.

###### Proof.

We show that the above $S^{W}$’s satisfy the Poisson equation (36) directly.

First, we consider the mean parameter $\mu_{i}$. Let $\bm{e}_{i}$ be the
$i$-th standard unit vector. From (58),

$\displaystyle\frac{\partial}{\partial\mu_{i}}\log
p(\bm{x},\bm{\theta})=-\frac{\partial}{\partial x_{i}}\log
p(\bm{x},\bm{\theta})=-\nabla_{\bm{x}}\log
p(\bm{x},\bm{\theta})^{\top}\bm{e}_{i}.$

Also, from Lemma 1,

$\displaystyle\nabla_{\bm{x}}(x_{i}-\mu_{i})=\bm{e}_{i},\quad\Delta_{\bm{x}}(x_{i}-\mu_{i})=0.$

Therefore,

$\displaystyle\nabla_{\bm{x}}\log
p({\bm{x}},{\bm{\theta}})^{\top}\nabla_{\bm{x}}(x_{i}-\mu_{i})+\Delta_{\bm{x}}(x_{i}-\mu_{i})+\frac{\partial}{\partial\mu_{i}}\log
p({\bm{x}},{\bm{\theta}})=0.$

Thus, the Wasserstein score function for the mean parameter $\mu_{i}$ is

$\displaystyle S_{\mu_{i}}^{{W}}(\bm{x},\bm{\theta})=x_{i}-\mu_{i}.$

Next, we consider the deformation parameter $\Lambda_{ij}$. Since

$\displaystyle\frac{\partial}{\partial\Lambda_{ij}}\|\Lambda(\bm{x}-\mu)\|$
$\displaystyle=\frac{1}{2}\|\Lambda(\bm{x}-\mu)\|^{-1}\frac{\partial}{\partial\Lambda_{ij}}(\bm{x}-\mu)^{\top}\Lambda^{2}(\bm{x}-\mu)$
$\displaystyle=\frac{1}{2}\|\Lambda(\bm{x}-\mu)\|^{-1}(\bm{x}-\mu)^{\top}\frac{\partial\Lambda^{2}}{\partial\Lambda_{ij}}(\bm{x}-\mu)$
$\displaystyle=\frac{1}{2}\|\Lambda(\bm{x}-\mu)\|^{-1}(\bm{x}-\mu)^{\top}(\Lambda\bm{e}_{i}\bm{e}_{j}^{\top}+\bm{e}_{i}\bm{e}_{j}^{\top}\Lambda)(\bm{x}-\mu),$

we have

$\displaystyle\frac{\partial}{\partial\Lambda_{ij}}\log p(\bm{x},\theta)$
$\displaystyle=\frac{\partial}{\partial\Lambda_{ij}}\log\det\Lambda+\frac{\partial}{\partial\Lambda_{ij}}\log
g(\|\Lambda(\bm{x}-\mu)\|)$
$\displaystyle=(\Lambda^{-1})_{ij}+\frac{g^{\prime}(\|\Lambda(\bm{x}-\mu)\|)}{g(\|\Lambda(\bm{x}-\mu)\|)}\frac{\partial}{\partial\Lambda_{ij}}\|\Lambda(\bm{x}-\mu)\|$
$\displaystyle=(\Lambda^{-1})_{ij}+\frac{g^{\prime}(\|\Lambda(\bm{x}-\mu)\|)}{2\|\Lambda(\bm{x}-\mu)\|g(\|\Lambda(\bm{x}-\mu)\|)}(\bm{x}-\mu)^{\top}(\Lambda\bm{e}_{i}\bm{e}_{j}^{\top}+\bm{e}_{i}\bm{e}_{j}^{\top}\Lambda)(\bm{x}-\mu)$
$\displaystyle=-\frac{g^{\prime}(\|\Lambda(\bm{x}-\mu)\|)}{\|\Lambda(\bm{x}-\mu)\|g(\|\Lambda(\bm{x}-\mu)\|)}\left(-\frac{1}{2}\bm{x}^{\top}L\bm{x}+\bm{x}^{\top}L\mu-\frac{1}{2}\mu^{\top}L\mu\right)+(\Lambda^{-1})_{ij},$

where
$L=\Lambda\bm{e}_{i}\bm{e}_{j}^{\top}+\bm{e}_{i}\bm{e}_{j}^{\top}\Lambda$. Let

$\displaystyle S(\bm{x})=\frac{1}{2}\bm{x}^{\top}A\bm{x}+b^{\top}\bm{x}-{\rm
E}_{\theta}\left[\frac{1}{2}\bm{x}^{\top}A\bm{x}+b^{\top}\bm{x}\right],$

where $A$ is the unique solution of the Sylvester equation
$\Lambda^{2}A+A\Lambda^{2}=-L$ and $b=-A\mu$. From Lemma 1 and Lemma 2,

$\displaystyle\Delta_{\bm{x}}S(\bm{x})={\rm tr}A=-\frac{1}{2}{\rm
tr}(\Lambda^{-2}L)=-\frac{1}{2}{\rm
tr}(\Lambda^{-2}(\Lambda\bm{e}_{i}\bm{e}_{j}^{\top}+\bm{e}_{i}\bm{e}_{i}^{\top}\Lambda))=-(\Lambda^{-1})_{ij}.$

Also,

$\displaystyle\nabla_{\bm{x}}\log
p({\bm{x}},{\bm{\theta}})^{\top}\nabla_{\bm{x}}S(\bm{x})$ $\displaystyle=$
$\displaystyle\frac{g^{\prime}(\|\Lambda(\bm{x}-\mu)\|)}{g(\|\Lambda(\bm{x}-\mu)\|)}\nabla_{\bm{x}}(\|\Lambda(\bm{x}-\mu)\|)^{\top}(A\bm{x}+b)$
$\displaystyle=$
$\displaystyle\frac{g^{\prime}(\|\Lambda(\bm{x}-\mu)\|)}{\|\Lambda(\bm{x}-\mu)\|g(\|\Lambda(\bm{x}-\mu)\|)}(\Lambda^{2}(\bm{x}-\mu))^{\top}(A\bm{x}+b)$
$\displaystyle=$
$\displaystyle\frac{g^{\prime}(\|\Lambda(\bm{x}-\mu)\|)}{\|\Lambda(\bm{x}-\mu)\|g(\|\Lambda(\bm{x}-\mu)\|)}\left(\frac{1}{2}\bm{x}^{\top}(\Lambda^{2}A+A\Lambda^{2})\bm{x}+\bm{x}^{\top}(\Lambda^{2}b-A\Lambda^{2}\mu)-\mu^{\top}\Lambda^{2}b\right)$
$\displaystyle=$
$\displaystyle\frac{g^{\prime}(\|\Lambda(\bm{x}-\mu)\|)}{\|\Lambda(\bm{x}-\mu)\|g(\|\Lambda(\bm{x}-\mu)\|)}\left(-\frac{1}{2}\bm{x}^{\top}L\bm{x}+\bm{x}^{\top}L\mu-\frac{1}{2}\mu^{\top}L\mu\right),$

where we used

$\displaystyle\nabla_{\bm{x}}(\|\Lambda(\bm{x}-\mu)\|)$
$\displaystyle=\frac{1}{2\|\Lambda(\bm{x}-\mu)\|}\nabla_{\bm{x}}(\|\Lambda(\bm{x}-\mu)\|^{2})$
$\displaystyle=\frac{1}{2\|\Lambda(\bm{x}-\mu)\|}\nabla_{\bm{x}}((\bm{x}-\mu)^{\top}\Lambda^{2}(\bm{x}-\mu))$
$\displaystyle=\frac{1}{\|\Lambda(\bm{x}-\mu)\|}\Lambda^{2}(\bm{x}-\mu).$

Therefore,

$\displaystyle\nabla_{\bm{x}}\log
p({\bm{x}},{\bm{\theta}})^{\top}\nabla_{\bm{x}}S(\bm{x})+\Delta_{\bm{x}}S(\bm{x})+\frac{\partial}{\partial\Lambda_{ij}}\log
p({\bm{x}},{\bm{\theta}})=0,$

which means that $S(\bm{x})$ is the Wasserstein score function for
$\Lambda_{ij}$.

∎

In Theorem 1, we considered $d^{2}$ elements of $\Lambda$ as independent
parameters. If we impose the symmetric constraint on $\Lambda$, then the
Wasserstein score function for the off-diagonal elements of $\Lambda$ becomes
the twice of the one in Theorem 1.

######  Corollary 1.

For the elliptically symmetric deformation model (58), the $W$-estimator
$\hat{\bm{\theta}}_{W}$ is the second-order moment estimator irrespective of
the waveform $f(\bm{z})=g(\|\bm{z}\|)$.

###### Proof.

From Theorem 1, the Wasserstein functions are quadratic functions of $\bm{x}$.
Thus, the estimating equation for the $W$-estimator is linear (non-singular)
with respect to the first and second-order empirical moments of $\bm{x}$.
Also, the number of estimated parameters is equal to the number of first and
second-order moments. Therefore, the $W$-estimator coincides with the second-
order moment estimator. ∎

Note that Gelbrich (1990) showed that the $L^{2}$-Wasserstein divergence for
the elliptically symmetric deformation model (58) does not depend on the
waveform and is given by

$D_{W}(p(\bm{x},\bm{\theta_{1}}),p(\bm{x},\bm{\theta_{2}}))=\|\mu_{1}-\mu_{2}\|^{2}+{\rm
tr}(\Lambda_{1}^{-2}+\Lambda_{2}^{-2}-2(\Lambda_{1}^{-1}\Lambda_{2}^{-2}\Lambda_{1}^{-1})^{1/2}).$

It is an interesting future problem to derive the Wasserstein score function
and $W$-estimator for general affine deformation models.

Regarding the geometric structure of the elliptically symmetric deformation
model (58), we obtain the following. See Figure 1.

######  Theorem 2.

When $f$ is spherically symmetric, the model $M_{f}$ is orthogonal to
${\cal{F}}_{S}$ at the origin ${\bm{\mu}}=0,\Lambda=I$ of $M_{f}$ with respect
to the Wasserstein metric.

###### Proof.

Let $\delta p({\bm{x}})$ be a tangent vector of ${\cal{F}}_{S}$ at the origin.
Since all $p({\bm{x}})$ in ${\cal{F}}_{S}$ satisfy the standard conditions
(1), (2), and (3), $\delta p({\bm{x}})$ is orthogonal to any quadratic
function of ${\bm{x}}$. Since the $W$-score functions of $M_{S}$ are quadratic
functions from Theorem 1, it implies that $\delta p({\bm{x}})$ is orthogonal
to the Waserstein functions of ${\bm{x}}$, which form the basis of the tangent
space of $M_{f}$, with respect to the Wasserstein metric. ∎

Note that $M_{f}$ is orthogonal to ${\cal{F}}_{S}$ only when $f$ is Gaussian
in the case of the Fisher-Rao metric.

Here, we discuss the relation between the current $W$-estimator and the
estimator (48). Amari and Matsuda (2022) studied the estimator
$\hat{{\bm{\theta}}}_{W}$ in (48) for the one-dimensional location-scale model
by using the order statistics $x_{(i)}$. This estimator minimizes the
Wasserstein distance between the empirical distribution and the model. Here,
we show that the estimator
$\hat{\bm{\theta}}_{W}=(\hat{{\bm{\mu}}}_{W},\hat{\sigma}_{W})$ given by (48)
is asymptotically equivalent to the $W$-estimator, which coincides with the
second-order moment estimator from Theorem 1. We assume $\mu=0$ without loss
of generality. The estimator (48) of the location is

$\hat{\mu}_{W}=\frac{1}{n}\sum_{i=1}^{n}x_{(i)}=\frac{1}{n}\sum_{i=1}^{n}x_{i},$
(61)

which is the empirical mean and coincides with the moment estimator. Also, the
estimator (48) of the scale is

$\hat{\sigma}_{W}=\sum_{i=1}^{n}k_{i}x_{(i)},$ (62)

where

$k_{i}=\int^{z_{i}}_{z_{i-1}}zf(z)dz.$ (63)

Here, $z_{i}$ is the $i$-th equipartition point of $f(z)$ defined by

$z_{i}=F^{-1}\left(\frac{i}{n}\right),$ (64)

where $F$ is the cumulative distribution function of $f(z)$. From $\mu=0$, we
have $x_{(i)}\approx\sigma z_{i}$ asymptotically. Hence,

$k_{i}\approx\frac{1}{n}z_{i}\approx\frac{1}{n}\frac{x_{(i)}}{\sigma},$ (65)

which leads to

$\hat{\sigma}_{W}=\sum_{i=1}^{n}k_{i}x_{(i)}\approx\frac{1}{n\sigma}\sum_{i=1}^{n}x^{2}_{(i)}=\frac{1}{n\sigma}\sum_{i=1}^{n}x^{2}_{i}.$
(66)

Since $\hat{\sigma}_{W}\approx\sigma$ asymptotically,

$\hat{\sigma}_{W}^{2}\approx\frac{1}{n}\sum_{i=1}^{n}x^{2}_{i}.$ (67)

This shows that
$\hat{\bm{\theta}}_{W}=(\hat{{\bm{\mu}}}_{W},\hat{\sigma}_{W})$ asymptotically
coincides with the second-order moment estimator.

##  6 W-efficiency implies robustness to waveform change

Following Li and Zhao (2023), we define the Wasserstein covariance
($W$-covariance) matrix ${\rm Var}_{\theta}^{\mathrm{W}}[\hat{\bm{\theta}}]$
of an estimator $\hat{\bm{\theta}}$ by the positive semidefinite matrix given
by

$\displaystyle{\rm Var}_{\theta}^{\mathrm{W}}[\hat{\bm{\theta}}]=({\rm
E}_{\theta}[(\nabla_{\bm{x}}\hat{\bm{\theta}}_{a})^{\top}(\nabla_{\bm{x}}\hat{\bm{\theta}}_{b})])_{ab}.$
(68)

Li and Zhao (2023) showed the Wasserstein–Cramer–Rao inequality

${\rm
Var}_{\theta}^{\mathrm{W}}(\hat{\bm{\theta}})\succeq\left(\frac{\partial}{\partial\theta}{\rm
E}_{\theta}[\hat{\bm{\theta}}]\right)^{\top}G_{\mathrm{W}}(\theta)^{-1}\left(\frac{\partial}{\partial\theta}{\rm
E}_{\theta}[\hat{\bm{\theta}}]\right),$ (69)

where

$\displaystyle\frac{\partial}{\partial\theta}{\rm
E}_{\theta}[\hat{\bm{\theta}}]:=\left(\frac{\partial}{\partial\theta_{j}}{\rm
E}_{\theta}[\hat{\bm{\theta}}_{i}]\right)_{ij}.$

A consistent estimator $\hat{\bm{\theta}}$ is said to be Wasserstein efficient
($W$-efficient) if its Wasserstein covariance asymptotically satisfies (69)
with equality. We give a proof of the Wasserstein–Cramer–Rao inequality based
on the Cauchy–Schwarz inequality in the Appendix.

We show that the $W$-covariance can be viewed as a measure of robustness of an
estimator to noise. Suppose that $X\sim p({\bm{x}},{\bm{\theta}})$ and we
estimate $\bm{\theta}$ from noisy observation $X+Z$ where $X$ and $Z$ are
independent, ${\rm E}[Z]=0$ and ${\rm Var}[Z]=\sigma^{2}I$ with $\sigma^{2}$
sufficiently small.

######  Theorem 3.

The Wasserstein covariance satisfies

$\displaystyle{\rm Var}^{\mathrm{W}}_{\theta}[\hat{\bm{\theta}}]_{ab}=$
$\displaystyle\lim_{\sigma^{2}\to 0}\frac{{\rm
Var}_{\theta}[\hat{\bm{\theta}}(X+Z)]_{ab}-{\rm
Var}_{\theta}[\hat{\bm{\theta}}(X)]_{ab}}{\sigma^{2}}$ (70)
$\displaystyle\quad-\frac{1}{2}\left({\rm
Cov}_{\theta}[\hat{\bm{\theta}}_{a}(X),\Delta\hat{\bm{\theta}}_{b}(X)]+{\rm
Cov}_{\theta}[\hat{\bm{\theta}}_{b}(X),\Delta\hat{\bm{\theta}}_{a}(X)]\right),$
(71)

where $\Delta$ is the Laplacian. In particular, when $\hat{\bm{\theta}}$ is
quadratic in $\bm{x}$,

${\rm Var}^{\mathrm{W}}_{\theta}[\hat{\bm{\theta}}]=\lim_{\sigma^{2}\to
0}\frac{{\rm Var}_{\theta}[\hat{\bm{\theta}}(X+Z)]-{\rm
Var}_{\theta}[\hat{\bm{\theta}}(X)]}{\sigma^{2}}.$ (72)

###### Proof.

By Taylor expansion, for sufficiently small $z$,

$\displaystyle\hat{\bm{\theta}}_{a}(x+z)\approx\hat{\bm{\theta}}_{a}(x)+\sum_{i}\frac{\partial\hat{\bm{\theta}}_{a}}{\partial
x_{i}}(x)z_{i}+\frac{1}{2}\sum_{i,j}\frac{\partial^{2}\hat{\bm{\theta}}_{a}}{\partial
x_{i}\partial x_{j}}(x)z_{i}z_{j}.$

From ${\rm E}[Z]=0$, ${\rm Var}[Z]=\sigma^{2}I$ and the independence of $X$
and $Z$,

$\displaystyle{\rm E}_{\theta}[\hat{\bm{\theta}}_{a}(X+Z)]$
$\displaystyle={\rm E}_{\theta}[\hat{\bm{\theta}}_{a}(X)]+\sum_{i}{\rm
E}_{\theta}\left[\frac{\partial\hat{\bm{\theta}}_{a}}{\partial
x_{i}}(X)\right]{\rm E}[z_{i}]+\frac{1}{2}\sum_{i,j}{\rm
E}_{\theta}\left[\frac{\partial^{2}\hat{\bm{\theta}}_{a}}{\partial
x_{i}\partial x_{j}}(X)\right]{\rm E}[z_{i}z_{j}]$ $\displaystyle={\rm
E}_{\theta}[\hat{\bm{\theta}}_{a}(X)]+\frac{1}{2}{\rm
E}_{\theta}[\Delta\hat{\bm{\theta}}_{a}(X)]\sigma^{2}.$

Also,

$\displaystyle{\rm
E}_{\theta}[\hat{\bm{\theta}}_{a}(X+Z)\hat{\bm{\theta}}_{b}(X+Z)]$
$\displaystyle={\rm
E}_{\theta}[\hat{\bm{\theta}}_{a}(X)\hat{\bm{\theta}}_{b}(X)]+\frac{1}{2}{\rm
E}_{\theta}[\hat{\bm{\theta}}_{a}(X)\Delta\hat{\bm{\theta}}_{b}(X)+\hat{\bm{\theta}}_{b}(X)\Delta\hat{\bm{\theta}}_{a}(X)]\sigma^{2}$
$\displaystyle\quad+{\rm
E}_{\theta}[(\nabla\hat{\bm{\theta}}_{a}(X))^{\top}(\nabla\hat{\bm{\theta}}_{b}(X))]\sigma^{2}+o(\sigma^{2})$
$\displaystyle={\rm
E}_{\theta}[\hat{\bm{\theta}}_{a}(X)\hat{\bm{\theta}}_{b}(X)]+\frac{1}{2}{\rm
E}_{\theta}[\hat{\bm{\theta}}_{a}(X)\Delta\hat{\bm{\theta}}_{b}(X)+\hat{\bm{\theta}}_{b}(X)\Delta\hat{\bm{\theta}}_{a}(X)]\sigma^{2}$
$\displaystyle\quad+{\rm
Var}_{\theta}[\hat{\bm{\theta}}]_{ab}\sigma^{2}+o(\sigma^{2}).$

Then,

$\displaystyle{\rm Var}_{\theta}[\hat{\bm{\theta}}(X+Z)]_{ab}$
$\displaystyle=$ $\displaystyle{\rm
E}_{\theta}[\hat{\bm{\theta}}_{a}(X+Z)\hat{\bm{\theta}}_{b}(X+Z)]-{\rm
E}_{\theta}[\hat{\bm{\theta}}_{a}(X+Z)]{\rm
E}_{\theta}[\hat{\bm{\theta}}_{b}(X+Z)]$ $\displaystyle=$ $\displaystyle{\rm
Var}_{\theta}[\hat{\bm{\theta}}(X)]_{ab}+{\rm
Var}^{\mathrm{W}}_{\theta}[\hat{\bm{\theta}}]_{ab}\sigma^{2}+\frac{1}{2}\left({\rm
Cov}_{\theta}[\hat{\bm{\theta}}_{a}(X),\Delta\hat{\bm{\theta}}_{b}(X)]+{\rm
Cov}_{\theta}[\hat{\bm{\theta}}_{b}(X),\Delta\hat{\bm{\theta}}_{a}(X)]\right)\sigma^{2}+o(\sigma^{2}),$

where the covariance term vanishes when $\hat{\bm{\theta}}$ is quadratic in
$\bm{x}$. ∎

Thus, the Wasserstein covariance quantifies the robustness of estimators
against the waveform change due to noise. In particular, from Corollary 1, the
Wasserstein covariance quantifies the robustness of the Wasserstein estimator
for elliptically symmetric deformation models. The Wasserstein–Cramer–Rao
inequality gives the limit of robustness. It is an interesting future problem
to investigate when the Wasserstein estimator attains the Wasserstein
efficiency. Note that the Fisher efficiency (in finite samples), which is
defined by the conventional Cramer–Rao inequality, is attained by the maximum
likelihood estimator if and only if the estimand is the expectation parameter
of an exponential family.

##  7 Contribution of waveform $f$ to $F$\- and $W$-efficiencies

We study how the waveform $f$ contributes to the $F$-efficiency and
$W$-efficiency of estimators. We first show the following theorem.

######  Theorem 4.

The Gaussian waveform $f$ has the following properties: When and only when $f$
is Gaussian, $F$-estimator $\hat{\bm{\theta}}_{F}$ and $W$-estimator
$\hat{\bm{\theta}}_{W}$ are identical, and they are both $F$\- and
$W$-efficient.

###### Proof.

For the standard Gaussian $f$,

$f({\bm{z}})=\frac{1}{(\sqrt{2x})^{d}}\exp\left\\{-\frac{|{\bm{z}}|^{2}}{2}\right\\},$
(73)

the $F$-score functions are

$S_{F}({\bm{x}},{\bm{\theta}})=\partial_{\bm{\theta}}\log
p({\bm{x}},{\bm{\theta}})=-{\bm{z}}\cdot\frac{\partial{\bm{z}}}{\partial{\bm{\theta}}}.$
(74)

Hence, the $F$-score functions are quadratic with respect to ${\bm{x}}$. So
they are equivalent to the $W$-score functions. On the contrary, when the
score functions are quadratic with respect to ${\bm{x}}$, the waveform $f$ is
Gaussian. ∎

When $f$ is not Gaussian, the $F$-efficiency of $\hat{\bm{\theta}}_{W}$ and
$W$-efficiency of $\hat{\bm{\theta}}_{F}$ degrade. When $f$ is close to
Gaussian, their cumulants of order larger than two are small. We use the Gram-
Charlie expansion of $f$ to represent how cumulants of the waveform $f$
contribute to the $F$-efficiency of $\hat{\bm{\theta}}_{W}$.

We study how waveform $f$ contributes to the amounts of Fisher information
$g_{F}$ and Wasserstein information $g_{W}$, when $f$ is close to the Gaussian
distribution. We use the Gram-Charlie expansion (McCullagh, 2018)

$f({\bm{x}})=S(\bm{x})\left\\{1+\frac{\kappa_{3}}{3!}\circ
h_{3}({\bm{x}})+\frac{\kappa_{4}}{4!}\circ h_{4}({\bm{x}})\right\\},$ (75)

where $\kappa_{i}$ are the $i$th order cumulant tensors and $h_{i}({\bm{x}})$
are the $i$th order tensorial Hermite polynomials, and $\circ$ denotes the
tensorial inner product such as

$\kappa_{3}\circ h_{3}=\sum\kappa_{3,ijk}h_{3}^{ijk}.$ (76)

Note that $\kappa_{1}=0,\kappa_{2}=1$ when $f$ is the standard Gaussian
distribution $\phi$,

$\phi({\bm{x}})=\frac{1}{(\sqrt{2\pi)}}\exp\left\\{-\frac{{\bm{x}}\cdot{\bm{x}}}{2}\right\\}.$
(77)

The logarithm of $p({\bm{x}},{\bm{\theta}})$ is expanded as

$l({\bm{x}},{\bm{\theta}})=-\log|\Lambda|+\left\\{-\frac{1}{2}|{\bm{z}}|^{2}+\frac{\kappa_{3}}{3!}\circ
h_{3}({\bm{z}})+\frac{\kappa_{4}}{4!}\circ h_{4}({\bm{z}})\right\\},$ (78)

in terms of ${\bm{z}}=\Lambda({\bm{x}}-{\bm{\mu}})$, where higher order terms
of $\kappa_{3}$ and $\kappa_{4}$ and terms of $\kappa_{i}$, $i\geq 5$ are
neglected.

Fisher information $g_{F}$ is given by

$g_{F}=-{\rm{E}}\left[\partial_{\bm{\theta}}\partial_{\bm{\theta}}l({\bm{x}},{\bm{\theta}})\right].$
(79)

We have

$\displaystyle\partial_{\bm{\theta}}l$ $\displaystyle=$
$\displaystyle\frac{\partial
l}{\partial{\bm{z}}}\cdot\frac{\partial{\bm{z}}}{\partial_{\bm{\theta}}}$ (80)
$\displaystyle\partial_{\bm{\theta}}\partial_{\bm{\theta}}l$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}l}{\partial{\bm{z}}\partial{\bm{z}}}\cdot\left(\frac{\partial{\bm{z}}}{\partial_{\bm{\theta}}}\frac{\partial{\bm{z}}}{\partial_{\bm{\theta}}}\right)+\frac{\partial
l}{\partial{\bm{z}}}\cdot\frac{\partial^{2}{\bm{z}}}{\partial_{\bm{\theta}}\partial_{\bm{\theta}}},$
(81)

where $\cdot$ denotes the inner product with respect to the indices of
${\bm{z}}$. From the derivative of Hermite polynomials, we have

$\displaystyle\partial_{\bm{z}}l$ $\displaystyle=$
$\displaystyle\Lambda^{-1}\left\\{-{\bm{z}}+\frac{\kappa_{3}}{2}\circ
h_{2}({\bm{z}})+\frac{\kappa_{4}}{6}\circ h_{3}({\bm{z}})\right\\}$ (82)
$\displaystyle\partial_{\bm{z}}\partial_{\bm{z}}l$ $\displaystyle=$
$\displaystyle\Lambda^{-1}\left\\{-I+\kappa_{3}\circ
h_{1}({\bm{z}})+\frac{\kappa_{4}}{2}\circ h_{2}({\bm{z}})\right\\},$ (83)

where $\kappa_{i}\circ h_{j}({\bm{z}})$ are tensorial polynomials of
${\bm{z}}$. On the other hand, $\partial{\bm{z}}/\partial{\bm{\theta}}$ are
given by

$\displaystyle\frac{\partial{\bm{z}}}{\partial{\bm{\mu}}}$ $\displaystyle=$
$\displaystyle\Lambda,$ (84)
$\displaystyle\frac{\partial{\bm{z}}}{\partial{\bm{\Lambda}}}$
$\displaystyle=$ $\displaystyle{\bm{x}}-{\bm{\mu}},$ (85)

In order to avoid complicated tensorial calculations, we study only the case
of $d=1$, that is the location-scale model. We show the results after simple
calculations.

$\displaystyle\frac{\partial^{2}l}{\partial\mu\partial\mu}$ $\displaystyle=$
$\displaystyle\lambda^{2}\left\\{-1+\kappa_{3}z+\frac{\kappa_{4}}{2}\left(z^{2}-1\right)\right\\}$
(86) $\displaystyle\frac{\partial^{2}l}{\partial\lambda\partial\mu}$
$\displaystyle=$
$\displaystyle\lambda^{2}\left\\{-2z+\frac{\kappa_{3}}{2}\left(3z^{2}-1\right)+\frac{\kappa_{4}}{3}\left(2z^{3}-3z\right)\right\\}$
(87) $\displaystyle\frac{\partial^{2}l}{\partial\lambda\partial\lambda}$
$\displaystyle=$
$\displaystyle\lambda^{2}\left\\{\left(-3z^{2}+1\right)+\frac{\kappa_{3}}{2}\left(4z^{3}-3z\right)+\frac{\kappa_{4}}{6}\left(4z^{4}-6z^{2}\right)\right\\}.$
(88)

Here, $z$ is subject to $f(z)$, so we have

$g_{F}(\mu,\lambda)=\lambda^{2}\left[\begin{array}[]{cc}1&-\kappa_{3}\\\
-\kappa_{3}&2-\kappa_{4}\end{array}\right].$ (89)

This shows how $g_{F}(\mu,\lambda)$ deviates from the Gaussian case depending
on $\kappa_{3}$ and $\kappa_{4}$.

It is also interesting to consider the case when $f$ has high-frequency wavy
structure. Since $F$-score functions are derivatives of the $\log$
probability, high-frequency components are sensitive to them, contributing to
the $F$-metric. However, by adding small noises to ${\bm{x}}$, those
components are smoothed out. Hence the $W$-metric is insensitive to the high-
frequency components.

We observe that, when $f$ includes a high frequency component such as

$f(x)=\phi({\bm{x}})\left\\{1+\varepsilon\sin\tau x\right\\},$ (90)

where $\varepsilon$ is small and $\tau$ is the frequency of small deviation,
$\partial^{2}l/\partial z^{2}$ has a component proportional to
$\varepsilon\tau^{2}$. Hence, the increment due to the high-frequency
component is proportional to $\tau^{2}$, implying that the increment of Fisher
information is proportional to $\tau^{2}$. Hence, high frequency ripples of
waveform $f$ increases $g_{F}$.

On the other hand, the $W$-estimator $\hat{\bm{\theta}}_{W}$ is not
$F$-efficient except for the Gaussian case. The loss of $F$-efficiency depends
on the waveform $f$. We again use the Gram-Charlie expansion and see the
effect of $\kappa_{3}$ and $\kappa_{4}$, assuming they are small. We show this
only in the location-scale model with $d=1$. Let us define the empirical
moments of order $s$ by

$m_{r}=\frac{1}{n}\sum^{n}_{i=1}x^{r}_{i},\quad r=1,2,\cdots$ (91)

$\hat{\bm{\theta}}_{W}$ uses only $m_{1}$ and $m_{2}$, discarding higher-order
moments. We calculate Fisher information $g_{F,W}({\bm{\theta}})$ included in
$\hat{\bm{\theta}}_{W}$. The Fisher information is covariance matrix of
$F$-score $\nabla l({\bm{x}},{\bm{\theta}})$. We classify
${\bm{x}}=\left(x_{1},\cdots,x_{n}\right)$ into classes specified by
$M=\left(m_{1},m_{2}\right)$, such that class $C_{M}$ consists of ${\bm{x}}$
having moments $s_{1},s_{2}$

$C_{M}=\left\\{{\bm{x}}\left|m_{1}({\bm{x}})=s_{1},m_{2}({\bm{x}})=s_{2}\right.\right\\}.$
(92)

Then, covariance of $\nabla l$ is decomposed into the sum of within-class
covariance and between class covariance,

${\rm{Cov}}[\nabla l]={\rm{E}}_{M}[{\rm{Cov}}[\nabla
l|M]]+{\rm{Cov}}[{\rm{E}}[\nabla l|M]],$ (93)

where ${\rm{E}}[\cdot|M]$ and ${\rm{Cov}}[\cdot|M]$ denotes the conditional
expectation and conditional covariance conditioned on $M$. Since
$\hat{\bm{\theta}}_{W}$ does not use higher-order moment information not
included in $M$, the Fisher information included in $\hat{\bm{\theta}}_{W}$ is
only the between-class covariance. The loss of information in
$\hat{\bm{\theta}}_{W}$ is

$\Delta g_{F,W}={\rm{E}}_{M}\left[{\rm{Cov}}[\nabla l|M]\right].$ (94)

The conditional expectation of $\nabla l$ is

$\displaystyle E\left[\partial_{\mu}l\left|s_{1},s_{2}\right.\right]$
$\displaystyle=$ $\displaystyle-
s_{1}+\frac{\kappa_{3}}{2}\left(s_{2}-1\right)+\frac{\kappa_{4}}{6}\left({\rm{E}}\left[s_{3}\left|M\right.\right]-3s_{1}\right)$
(95) $\displaystyle E\left[\partial_{\lambda}l\left|s_{1},s_{2}\right.\right]$
$\displaystyle=$ $\displaystyle-
s_{2}+\frac{\kappa_{3}}{2}\left\\{\left({\rm{E}}\left[s_{3}\left|M\right.\right]\right)-s_{1}\right\\}+\frac{\kappa_{4}}{6}\left\\{{\rm{E}}\left[s_{4}\left|M\right.\right]-3s_{2}\right\\}.$
(96)

Hence, the conditional covariances are

$\displaystyle{\rm{Cov}}\left[\partial_{\mu}l\left|M\right.\right]$
$\displaystyle=$
$\displaystyle\frac{\kappa^{2}_{4}}{36}{\rm{Cov}}\left[s_{3}\left|M\right.\right],$
(97) $\displaystyle{\rm{Cov}}\left[\partial_{\lambda}l\left|M\right.\right]$
$\displaystyle=$
$\displaystyle\frac{\kappa^{2}_{3}}{4}{\rm{Cov}}\left[s_{3}\left|M\right.\right]+\frac{\kappa^{2}_{4}}{36}{\rm{Cov}}\left[s_{4}\left|M\right.\right].$
(98)

${\rm{Cov}}\left[\partial_{\mu}l,\partial_{\lambda}l\left|M\right.\right]=\frac{\kappa_{3}\kappa_{4}}{12}{\rm{E}}\left[s_{3}|M\right]{\rm{E}}\left[s_{4}|M\right].$
(99)

It should be noted that $s_{3}$ and $s_{4}$ are asymptotically independent of
$M$, because $\left(s_{1},s_{2},s_{3},s_{4}\right)$ are jointly Gaussian and
asymtptically independent.

We thus have asymptotically

$\displaystyle{\rm{E}}\left[s_{3}\left|M\right.\right]$ $\displaystyle=$
$\displaystyle{\rm{E}}\left[x^{3}\right]=\frac{\kappa_{3}}{\lambda^{3}},$
(100) $\displaystyle{\rm{E}}\left[s_{4}\left|M\right.\right]$ $\displaystyle=$
$\displaystyle{\rm{E}}\left[x^{4}\right]=\frac{3\left(1+\kappa_{4}\right)}{\lambda^{4}}+\frac{6}{\lambda^{2}}\mu^{2}+\mu^{4},$
(101)

and so on. Summing up all these results, we have $\Delta g_{F,W}$ in terms of
$\kappa_{3}$ and $\kappa_{4}$.

##  8 Conclusion

Statistical inference has so far been studied mostly based on information
geometry from the Fisherian point of view, with remarkable success. It is
based on the likelihood principle, and the invariant divergence has played a
fundamental role. However, Wasserstein divergence gives another viewpoint,
which is based on the geometric structure of the sample space $X$. There are
many applications of the Wasserstein geometry not only to the transportation
problem but to vision analysis, signal analysis and AI in which the geometry
of $X$ is sensible.

We studied the Wasserstein statistics using the framework of Li and Zhao
(2023), proving that the Wasserstein covariance quantifies robustness against
the convolutional waveform deformation due to observation noise. We further
studied $W$-statistics of the affine deformation model. We showed
$F$-efficiency and $W$-efficiency of estimators $\hat{\bm{\theta}}_{F}$ and
$\hat{\bm{\theta}}_{W}$. We elucidated how the waveform $f$ contributes to the
efficiencies. The Gaussian distribution gives the only waveform in which the
$F$-estimator and $W$-estimator coincide, satisfying both efficiencies.

The present paper is only a first step to construct general Wasserstein
statistics. In future work, we need to use more general statistical models. We
also need to extend our approach to statistical theories of hypothesis
testing, pattern classification, clustering and many other statistical
problems based on the Wasserstein geometry.

## Acknowledgements

We thank Asuka Takatsu and Tomonari Sei for helpful comments. We thank Emi
Namioka for drawing the figures. Takeru Matsuda was supported by JSPS KAKENHI
Grant Numbers 19K20220, 21H05205, 22K17865 and JST Moonshot Grant Number
JPMJMS2024.

## References

  * Amari (2016) Amari, S. (2016). Information Geometry and Its Applications. Springer.
  * Amari et al. (2018) Amari, S., Karakida, R. & Oizumi, M. (2018). Information geometry connecting Wasserstein distance and Kullback–Leibler divergence via the entropy-relaxed transportation problem. Information Geometry, 1, 13–37.
  * Amari et al. (2019) Amari, S., Karakida, R., Oizumi, M. & Cuturi, M. (2019). Information geometry for regularized optimal transport and barycenters of patterns. Neural Computation, 31, 827–848.
  * Amari and Matsuda (2022) Amari, S. & Matsuda, T. (2022). Wasserstein statistics in one-dimensional location scale models Annals of the Institute of Statistical Mathematics, 74, 33–47.
  * Amari and Nagaoka (2007) Amari, S. & Nagaoka, H. (2016). Methods of Information Geometry. American Mathematical Society.
  * Arjovsky et al. (2017) Arjovsky, M., Chintala, S. & Bottou, L. (2017). Wasserstein GAN. arXiv:1701.07875.
  * Ay et al. (2017) Ay, N. and Jost, J. and Vân Lê, H. & Schwachhöfer, L. (2017). Information Geometry. Springer.
  * Bassetti et al. (2006) Bassetti, F., Bodini, A. & Regazzini, E. (2006). On minimum Kantorovich distance estimators. Statistics & Probability Letters, 76, 1298–1302.
  * Benamou and Brenier (2000) Benamou, J. D. & Brenier, Y. (2000). A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numerische Mathematik, 84, 375–393.
  * Bernton et al. (2019) Bernton, E., Jacob, P. E., Gerber, M. & Robert, C. P. (2019). On parameter estimation with the Wasserstein distance. Information and Inference: A Journal of the IMA, 8, 657–676.
  * Bhatia (1997) Bhatia, R. (1997). Matrix Analysis. Springer.
  * Brenier (1999) Brenier, Y. (1999). Minimal geodesics on groups of volume-preserving maps and generalized solutions of the Euler equations. Comm. Pure Appl. Math., 52, 411–452.
  * Chentsov (1982) Chentsov, N. (1982). Statistical decision rules and optimal inference. American Mathematical Society.
  * Chen et al. (2021) Chen, Y., Lin, Z. & Müller, H. G. (2021). Wasserstein regression. Journal of the American Statistical Association, accepted.
  * Fronger et al. (2015) Fronger, C., Zhang, C., Mobahi, H., Araya-Polo, M. & Poggio, T. (2015). Learning with a Wasserstein loss. Advances in Neural Information Processing Systems 28 (NIPS 2015).
  * Gelbrich (1990) Gelbrich, M. (1990). On a formula for the L2 Wasserstein metric between measures on Euclidean and Hilbert spaces. Mathematics Nachrichten 147, 185–203.
  * Givens and Shortt (1984) Givens, C. R. & Shortt, R. M. (1984). A class of Wasserstein metrics for probability distributions. Michigan Mathematics Journal 31, 231–240.
  * Imaizumi et al. (2022) Imaizumi, M., Ota, H. & Hamaguchi, T. (2022). Hypothesis test and confidence analysis with Wasserstein distance on general dimension. Neural Computation 34. 1448–1487.
  * Ito (2023) Ito, S. (2023). Geometric thermodynamics for the Fokker–Planck equation: stochastic thermodynamic links between information geometry and optimal transport. Information Geometry.
  * Khan et al. (2022) Khan, G. & Zhang, J. (2022). When optimal transport meats information geometry. Infromation Geometry, 1–32.
  * Kurose et al. (2019) Kurose, T., Yoshizawa, S. & Amari, S. (2019). Optimal transportation plan with generalized entropy regularization. submitted.
  * Li and Montufar (2020) Li, W. & Montúfar, G. (2020). Ricci curvature for parametric statistics via optimal transport. Information Geometry, 3, 89-–117.
  * Li and Zhao (2023) Li, W. & Zhao, J. (2023). Wasserstein information matrix. Information Geometry, accepted.
  * Matsuda and Strawderman (2021) Matsuda, T. & Strawderman, W. E. (2021). Predictive density estimation under the Wasserstein loss. Journal of Statistical Planning and Inference, 210, 53–63.
  * McCullagh (2018) McCullagh, P. (2018). Tensor Methods in Statistics. Dover.
  * Montavon et al. (2015) Montavon, G., Müller, K. R. & Cuturi, M. (2015). Wasserstein training for Boltzmann machine. Advances in Neural Information Processing Systems 29 (NIPS 2016).
  * Ollila and Tyler (2014) Ollila, E. & Tyler, D. (2014). Regularized $M$-estimators of scatter matrix. IEEE Transactions on Signal Processing, 62, 6059–6070.
  * Otto (2001) Otto, F. (2001). The geometry of dissipative evolution equations: The porous medium equation. Commun. Partial Differ. Equations, 26, 101–174.
  * Panaretos and Zemel (2019) Panaretos, V. M. & Zemel, Y. (2019). Statistical aspects of Wasserstein distances. arXiv:1806.05500v3.
  * Panaretos and Zemel (2022) Panaretos, V. M. & Zemel, Y. (2022). An Invitation to Statistics in Wasserstein Space. Springer.
  * Peyré and Cuturi (2019) Peyré, G. & Cuturi, M. (2019). Computational optimal transport: With Applications to Data Science. Foundations and Trends® in Machine Learning, 11, 355–607.
  * Rankin and Wong (2023) Rankin, C. & Wong, TK, L. (2023). Bregman-Wasserstein divergence: geometry and applications. arXiv:2302.05833.
  * Santambrogio (2015) Santambrogio, F. (2015). Optimal Transport for Applied Mathematicians. Springer.
  * Takatsu (2011) Takatsu, A. (2011). Wasserstein geometry of Gaussian measures. Osaka Journal of Mathematics, 48, 1005–1026.
  * Takatsu and Yokota (2012) Takatsu, A. & Yokota, T.(2012). Cone structure of $L^{2}$-Wasserstein spaces. Journal of Topology and Analysis, 4, 237–253.
  * Tyler (1987) Tyler, D. (1987). A distribution-free M-estimator of multivariate scatter. Annals of Statistics, 15, 234–251.
  * van der Vaart (1998) van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press.
  * Villani (2003) Villani, C. (2003). Topics in Optimal Transportation. American Mathematical Society.
  * Villani (2009) Villani, C. (2009). Optimal Transport: Old and New. Springer.
  * Wang and Li (2020) Wang, Y. & Li, W. (2020). Information Newton’s flow: Second-order optimization method in probability space. arXiv:2001.04341.
  * Yatracos (2022) Yatracos, Y. G. (2022). Limitations of the Wasserstein MDE for univariate data. Statistics and Computing volume, 32, 32–95.

##  Appendix A Proof of Wasserstein–Cramer–Rao inequality

For random vectors $A$ and $B$,

$\displaystyle 0\leq{\rm E}\|tA+B\|^{2}={\rm E}[\|A\|^{2}]t^{2}+2{\rm
E}[A^{\top}B]t+{\rm E}[\|B\|^{2}]$ (102)

for every $t$. Thus, by considering the discriminant of the quadratic
equation,

$\displaystyle{\rm E}[A^{\top}B]^{2}\leq{\rm E}[\|A\|^{2}]{\rm E}[\|B\|^{2}].$
(103)

Substituting

$\displaystyle A=\sum_{i}a_{i}\nabla_{\bm{x}}\hat{\theta}_{i},\quad
B=\sum_{j}b_{j}\nabla_{\bm{x}}S_{j}^{W}$ (104)

into (103) yields

$\displaystyle\left(\sum_{i,j}a_{i}b_{j}{\rm
E}_{\theta}[(\nabla_{\bm{x}}\hat{\theta}_{i})^{\top}(\nabla_{\bm{x}}S_{j}^{W})]\right)^{2}\leq{\rm
E}_{\theta}\left[\left\|\sum_{i}a_{i}\nabla_{\bm{x}}\hat{\theta}_{i}\right\|^{2}\right]{\rm
E}_{\theta}\left[\left\|\sum_{j}b_{j}\nabla_{\bm{x}}S_{j}^{W}\right\|^{2}\right].$
(105)

From the property (36) of the Wasserstein score function,

$\displaystyle\frac{\partial}{\partial\theta_{j}}{\rm
E}_{\theta}[\hat{\theta}_{i}]$
$\displaystyle=\int\hat{\theta}_{i}(x)\frac{\partial}{\partial\theta_{j}}p(x,\theta){\rm
d}x$ (106)
$\displaystyle=-\int\hat{\theta}_{i}(x)\nabla_{\bm{x}}\cdot(p(x,\theta)\nabla_{\bm{x}}S_{j}^{W}(x,\theta)){\rm
d}x$ (107)
$\displaystyle=-\int(\nabla_{\bm{x}}\cdot(\hat{\theta}_{i}(x)p(x,\theta)\nabla_{\bm{x}}S_{j}^{W}(x,\theta))-(\nabla_{\bm{x}}\hat{\theta}_{i}(x))^{\top}(\nabla_{\bm{x}}S_{j}^{W}(x,\theta))p(x,\theta)){\rm
d}x$ (108) $\displaystyle={\rm
E}_{\theta}[(\nabla_{\bm{x}}\hat{\theta}_{i})^{\top}(\nabla_{\bm{x}}S_{j}^{W})],$
(109)

which yields

$\displaystyle\sum_{i,j}a_{i}b_{j}{\rm
E}_{\theta}[(\nabla_{\bm{x}}\hat{\theta}_{i})^{\top}(\nabla_{\bm{x}}S_{j}^{W})]=\sum_{i,j}a_{i}b_{j}\frac{\partial}{\partial\theta_{j}}{\rm
E}_{\theta}[\hat{\theta}_{i}]=a^{\top}\left(\frac{\partial}{\partial\theta}{\rm
E}_{\theta}[\hat{\theta}]\right)b.$ (110)

From the definition of the Wasserstein covariance (68),

$\displaystyle{\rm
E}_{\theta}\left[\left\|\sum_{i}a_{i}\nabla_{\bm{x}}\hat{\theta}_{i}\right\|^{2}\right]$
$\displaystyle={\rm
E}_{\theta}\left[\sum_{i,j}a_{i}a_{j}(\nabla_{\bm{x}}\hat{\theta}_{i})^{\top}(\nabla_{\bm{x}}\hat{\theta}_{j})\right]$
(111) $\displaystyle=\sum_{i,j}a_{i}a_{j}{\rm
Var}_{\theta}^{\mathrm{W}}(\hat{\theta})_{ij}$ (112)
$\displaystyle=a^{\top}{\rm Var}_{\theta}^{\mathrm{W}}(\hat{\theta})a.$ (113)

From the definition of the Wasserstein information matrix (42),

$\displaystyle{\rm
E}_{\theta}\left[\left\|\sum_{j}b_{j}\nabla_{\bm{x}}S_{j}^{W}\right\|^{2}\right]$
$\displaystyle={\rm
E}_{\theta}\left[\sum_{i,j}b_{i}b_{j}(\nabla_{\bm{x}}S_{i}^{W})^{\top}(\nabla_{\bm{x}}S_{j}^{W})\right]$
(114) $\displaystyle=\sum_{i,j}b_{i}b_{j}G_{\mathrm{W}}(\theta)_{ij}$ (115)
$\displaystyle=b^{\top}G_{\mathrm{W}}(\theta)b.$ (116)

Substituting (110), (113) and (116) into (105), we obtain

$\displaystyle\left(a^{\top}\left(\frac{\partial}{\partial\theta}{\rm
E}_{\theta}[\hat{\theta}]\right)b\right)^{2}\leq a^{\top}{\rm
Var}_{\theta}^{\mathrm{W}}(\hat{\theta})a\cdot
b^{\top}G_{\mathrm{W}}(\theta)b.$ (117)

By putting

$\displaystyle
b=G_{\mathrm{W}}(\theta)^{-1}\left(\frac{\partial}{\partial\theta}{\rm
E}_{\theta}[\hat{\theta}]\right)a,$ (118)

it leads to

$\displaystyle a^{\top}\left(\frac{\partial}{\partial\theta}{\rm
E}_{\theta}[\hat{\theta}]\right)G_{\mathrm{W}}(\theta)^{-1}\left(\frac{\partial}{\partial\theta}{\rm
E}_{\theta}[\hat{\theta}]\right)a\leq a^{\top}{\rm
Var}_{\theta}^{\mathrm{W}}(\hat{\theta})a,$ (119)

which is equal to the Wasserstein Cramer–Rao inequality (69).