# Entropy-regularized optimal transport on multivariate normal and $q$-normal
distributions

QIJUN<EMAIL_ADDRESS>[ KEI
<EMAIL_ADDRESS>[ Department of Mathematics, Faculty of
Science and Technology, Keio University, Yokohama, Japan.  Department of
Mathematics, Faculty of Science and Technology, Keio University, Yokohama,
Japan.

###### Abstract

The distance and divergence of the probability measures play a central role in
statistics, machine learning, and many other related fields. The Wasserstein
distance has received much attention in recent years because of its
distinctions from other distances or divergences. Although computing the
Wasserstein distance is costly, entropy-regularized optimal transport was
proposed to computationally efficiently approximate the Wasserstein distance.
The purpose of this study is to understand the theoretical aspect of entropy-
regularized optimal transport. In this paper, we focus on entropy-regularized
optimal transport on multivariate normal distributions and $q$-normal
distributions. We obtain the explicit form of the entropy-regularized optimal
transport cost on multivariate normal and $q$-normal distributions; this
provides a perspective to understand the effect of entropy regularization,
which was previously known only experimentally. Furthermore, we obtain the
entropy-regularized Kantorovich estimator for the probability measure that
satisfies certain conditions. We also demonstrate how the Wasserstein
distance, optimal coupling, geometric structure, and statistical efficiency
are affected by entropy regularization in some experiments. In particular, our
results about the explicit form of the optimal coupling of the Tsallis
entropy-regularized optimal transport on multivariate $q$-normal distributions
and the entropy-regularized Kantorovich estimator are novel and will become
the first step towards the understanding of a more general setting.

Optimal trainsport,

Wasserstein distance,

Entropy regularization,

Tsallis entropy,

$q$-normal distribution,

###### keywords:

and

##  1 Introduction

Comparing probability measures is a fundamental problem in statistics and
machine learning. A classical way to compare probability measures is the
Kullback–Leibler divergence. Let $M$ be a measurable space and $\mu,\nu$ be
the probability measure on $M$; then, the Kullback–Leibler divergence is
defined as:

$\mathrm{KL}(\mu|\nu)=\int_{M}d\mu\log\frac{d\mu}{d\nu}.$ (1)

The Wasserstein distance [53], also known as the earth mover distance [47], is
another way of comparing probability measures. It is a metric on the space of
probability measures derived by the mass transportation theory of two
probability measures. Informally, optimal transport theory considers an
optimal transport plan between two probability measures under a cost function,
and the Wasserstein distance is defined by the minimum total transport cost. A
significant difference between the Wasserstein distance and the
Kullback–Leibler divergence is that the former can reflect the metric
structure, whereas the latter cannot. The Wasserstein distance can be written
as:

$W_{p}(\mu,\nu):=\left\\{\inf_{\pi\in\Pi(\mu,\nu)}\int_{M\times
M}d(x,y)^{p}d\pi(x,y)\right\\}^{\frac{1}{p}},$ (2)

where $d(\cdot,\cdot)$ is a distance function on a measurable metric space $M$
and $\Pi(\mu,\nu)$ denotes the set of probability measures on $M\times M$,
whose marginal measures correspond to $\mu$ and $\nu$. In recent years, the
application of optimal transport and the Wasserstein distance has been studied
in many fields such as statistics, machine learning, and image processing. For
example, Reference [49] generated the interpolation of various three-
dimensional (3D) objects using the Wasserstein barycenter. In the field of
word embedding in natural language processing, Reference [40] embedded each
word as an elliptical distribution, and the Wasserstein distance was applied
between the elliptical distributions. There are many studies on the
applications of optimal transport to deep learning, including [21, 6, 42].
Moreover, Reference [50] analyzed the denoising autoencoder [54] with gradient
flow in the Wasserstein space.

In the application of the Wasserstein distance, it is often considered in a
discrete setting where $\mu$ and $\nu$ are discrete probability measures.
Then, obtaining the Wasserstein distance between $\mu$ and $\nu$ can be
formulated as a linear programming problem. In general, however, it is
computationally intensive to solve such linear problems and obtain the optimal
coupling of two probability measures. For such a situation, a novel numerical
method, entropy regularization, was proposed by [16],

$C_{\lambda}(\mu,\nu):=\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}c(x,y)\pi({x},{y})d{x}d{y}-\lambda\mathrm{Ent}(\pi).$
(3)

This is a relaxed formulation of the original optimal transport of a cost
function $c(\cdot,\cdot)$, in which the negative Shannon entropy
$-\mathrm{Ent}(\cdot)$ is used as a regularizer. For a small $\lambda$,
$C_{\lambda}(\mu,\nu)$ can approximate the $p$-th power of the Wasserstein
distance between two discrete probability measures, and it can be computed
efficiently by using Sinkhorn’s algorithm [48].

More recently, many studies have been published on improving the computational
efficiency. According to [33], the most computationally efficient algorithm at
this moment to solve the linear problem for the Wasserstein distance is
Lee–Sidford linear solver [30], which runs in $O(n^{2.5})$. Reference [19]
proved that a complexity bound for the Sinkhorn algorithm is
$\tilde{O}(n^{2}\varepsilon^{-2})$, where $\varepsilon$ is the desired
absolute performance guarantee. After [16] appeared, various algorithms have
been proposed. For example, Reference [7] adopted stochastic optimization
schemes for solving the optimal transport. The Greenkhorn algorithm [3] is the
greedy variant of the Sinkhorn algorithm, and Reference [33] proposed its
acceleration. Many other approaches such as adapting a variety of standard
optimization algorithms to approximate the optimal transport problem can be
found in [11, 17, 33, 32]. Several specialized Newton-type algorithms [2, 13]
achieve complexity bound $\tilde{O}(n^{2}\varepsilon^{-1})$ [10, 45], which
are the best ones in terms of computational complexity at the present moment.

Moreover, entropy-regularized optimal transport has another advantage. Because
of the differentiability of the entropy-regularized optimal transport and the
simple structure of Sinkhorn’s algorithm, we can easily compute the gradient
of the entropy-regularized optimal transport cost and optimize the parameter
of a parametrized probability distribution by using numerical differentiation
or automatic differentiation. Then, we can define a differentiable loss
function that can be applied to various supervised learning methods [20].
Entropy-regularized optimal transport can be used to approximate not only the
Wasserstein distance, but also its optimal coupling as a mapping function.
Reference [15] adopted the optimal coupling of the entropy-regularized optimal
transport as a mapping function from one domain to another.

Despite the empirical success of the entropy-regularized optimal transport,
its theoretical aspect is less understood. Reference [31] studied the expected
Wasserstein distance between a probability measure and its empirical version.
Similarly, Reference [38] showed the consistency of the entropy-regularized
optimal transport cost between two empirical distributions. Reference [46]
showed that minimizing the entropy-regularized optimal transport cost between
empirical distributions is equivalent to a type of maximum likelihood
estimator. Reference [8] considered Wasserstein generative adversarial
networks with an entropy regularization. Reference [4] constructed information
geometry from the convexity of the entropy-regularized optimal transport cost.

Our intrinsic motivation of this study is to produce an analytical solution
about the entropy-regularized optimal transport problem between continuous
probability measures so that we can gain insight into the effects of entropy
regularization in a theoretical, as well as an experimental way. In our study,
we generalized the Wasserstein distance between two multivariate normal
distributions by entropy regularization. We derived the explicit form of the
entropy-regularized optimal transport cost and its optimal coupling, which can
be used to analyze the effect of entropy regularization directly. In general,
the nonregularized Wasserstein distance between two probability measures and
its optimal coupling cannot be expressed in a closed form; however, Reference
[18] proved the explicit formula for multivariate normal distributions.
Theorem 3.3 is a generalized form of [18]. We obtain an explicit form of the
entropy-regularized optimal transport between two multivariate normal
distributions. Furthermore, by adopting the Tsallis entropy [52] as the
entropy regularization instead of the Shannon entropy, our theorem can be
generalized to multivariate $q$-normal distributions.

Some readers may find it strange to study the entropy-regularized optimal
transport for multivariate normal distributions, where the exact
(nonregularized) optimal transport has been obtained explicitly. However, we
think it is worth studying from several perspectives:

  * •

Normal distributions are the simplest and best-studied probability
distributions, and thus, it is useful to examine the regularization
theoretically in order to infer results for other distributions. In
particular, we will partly answer the questions “How much do entropy
constraints affect the results?” and “What does it mean to constrain by the
entropy?” for the simplest cases. Furthermore, as a first step in constructing
a theory for more general probability distributions, in Section 4, we propose
a generalization to multivariate $q$-normal distributions.

  * •

Because normal distributions are the limit distributions in asymptotic
theories using the central limit theorem, studying normal distributions is
necessary for the asymptotic theory of regularized Wasserstein distances and
estimators computed by them. Moreover, it was proposed to use the entropy-
regularized Wasserstein distance to compute a lower bound of the
generalization error for a variational autoencoder [8]. The study of the
asymptotic behavior of such bounds is one of the expected applications of our
results.

  * •

Though this has not yet been proven theoretically, we suspect that entropy
regularization is efficient not only for computational reasons, such as the
use of the Sinkhorn algorithm, but also in the sense of efficiency in
statistical inference. Such a phenomenon can be found in some existing
studies, including [5]. Such statistical efficiency is confirmed by some
experiments in Section 6.

The remainder of this paper is organized as follows. First, we review some
definitions of optimal transport and entropy regularization in Section 2.
Then, in Section 3, we provide an explicit form of the entropy-regularized
optimal transport cost and its optimal coupling between two multivariate
normal distributions. We also extend this result to $q$-normal distributions
for Tsallis entropy regularization in Section 4. In Section 5, we obtain the
entropy-regularized Kantorovich estimator of probability measures on
$\mathbb{R}^{n}$ with a finite second moment that are absolutely continuous
with respect to the Lebesgue measure in Theorem 5.1. We emphasize that Theorem
5.1 is not limited to the case of multivariate normal distribution, but can
handle a wider range of probability measures. We analyze how entropy
regularization affects the optimal result experimentally in certain sections.

We note that after publishing the preprint version of the paper, we found
closely related results [24, 34] reported within half a year. In Janati et al.
[24], they proved the same result as Theorem 3.3 based on solving the fixed-
point equation behind Sinkhorn’s algorithm. Their results include the
unbalanced optimal transport between unbalanced multivariate normal
distributions. They also studied the convexity and differentiability of the
objective function of the entropy-regularized optimal transport. In [34], the
same closed-form as Theorem 3.3 was proven by ingeniously using the
Schrödinger system. Although there are some overlaps, our paper has
significant novelty in the following respects. Our proof is more direct than
theirs and can be extended directly to the proof for the Tsallis entropy-
regularized optimal transport between multivariate $q$-normal distributions
provided in Section 4. Furthermore, Corollaries 3.3.1 and 3.3.2 are novel and
important results to evaluate how much the entropy regularization affects the
estimation results or not at all. We also obtain the entropy-regularized
Kantorovich estimator in Theorem 5.1.

##  2 Preliminary

In this section, we review some definitions of optimal transport and entropy-
regularized optimal transport. These definitions were referred to in [44, 53].
In this section, we use a tuple $(M,\Sigma)$ for a set $M$ and
$\sigma$-algebra on $M$ and $\mathcal{P}(X)$ for the set of all probability
measures on a measurable space $X$.

######  Definition 2.1 (Pushforward measure).

Given measurable spaces $(M_{1},\Sigma_{1})$ and $(M_{2},\Sigma_{2})$, a
measure $\mu:\Sigma_{1}\rightarrow[0,+\infty]$, and a measurable mapping
$\varphi:M_{1}\rightarrow M_{2}$, the pushforward measure of $\mu$ by
$\varphi$ is defined by:

$\forall B\in\Sigma_{2},\
\varphi_{\\#}\mu(B):=\mu\left(\varphi^{-1}(B)\right).$ (4)

######  Definition 2.2 (Optimal transport map).

Consider a measurable space $(M,\Sigma)$, and let $c:M\times
M\rightarrow\mathbb{R}_{+}$ denote a cost function. Given
$\mu,\nu\in\mathcal{P}(M)$, we call $\varphi:M\rightarrow M$ the optimal
transport map if $\varphi$ realizes the infimum of:

$\inf_{\varphi_{\\#}\mu=\nu}\int_{M}c(x,\varphi(x))d\mu(x).$ (5)

This problem was originally formalized by [39]. However, the optimal transport
map does not always exist. Then, Kantorovich considered a relaxation of this
problem in [27].

######  Definition 2.3 (Coupling).

Given $\mu,\nu\in\mathcal{P}(M)$, the coupling of $\mu$ and $\nu$ is a
probability measure on $M\times M$ that satisfies:

$\forall A\in\Sigma,\ \pi(A\times M)=\mu(A),\quad\pi(M\times A)=\nu(A).$ (6)

######  Definition 2.4 (Kantorovich problem).

The Kantorovich problem is defined as finding a coupling $\pi$ of $\mu$ and
$\nu$ that realizes the infimum of:

$\int_{M\times M}c(x,y)d\pi(x,y).$ (7)

Hereafter, let $\Pi(\mu,\nu)$ be the set of all couplings of $\mu$ and $\nu$.
When we adopt a distance function as the cost function, we can define the
Wasserstein distance.

######  Definition 2.5 (Wasserstein distance).

Given $p\geq 1$, a measurable metric space $(M,\Sigma,d)$, and
$\mu,\nu\in\mathcal{P}(M)$ with a finite $p$-th moment, the $p$-Wasserstein
distance between $\mu$ and $\nu$ is defined as:

$W_{p}(\mu,\nu):=\inf_{\pi\in\Pi(\mu,\nu)}\left(\int_{M\times
M}d(x,y)^{p}d\pi(x,y)\right)^{\frac{1}{p}}.$ (8)

Now, we review the definition of entropy-regularized optimal transport on
$\mathbb{R}^{n}$.

######  Definition 2.6 (Entropy-regularized optimal transport).

Let $\mu,\nu\in\mathcal{P}(\mathbb{R}^{n})$, $\lambda>0$, and let $\pi(x,y)$
be the density function of the coupling of $\mu$ and $\nu$, whose reference
measure is the Lebesgue measure. We define the entropy-regularized optimal
transport cost as:

$C_{\lambda}(\mu,\nu):=\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}c(x,y)\pi({x},{y})d{x}d{y}-\lambda\mathrm{Ent}(\pi),$
(9)

where $\mathrm{Ent}(\cdot)$ denotes the Shannon entropy of a probability
measure:

$\mathrm{Ent}(\pi)=-\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}\pi(x,y)\log\pi(x,y)dxdy.$
(10)

There is another variation in entropy-regularized optimal transport defined by
the relative entropy instead of the Shannon entropy:

$\tilde{C}_{\lambda}(\mu,\nu):=\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}c(x,y)\pi({x},{y})d{x}d{y}+\lambda\mathrm{KL}(\pi|d\mu\otimes
d\nu).$ (11)

This is definable even when $\Pi(\mu,\nu)$ includes a coupling that is not
absolutely continuous with respect to the Lebesgue measure. We note that when
both $\mu$ and $\nu$ are absolutely continuous, the infimum is attained by the
same $\pi$ for $C_{\lambda}$ and $\tilde{C}_{\lambda}$, and it depends only on
$\mu$ and $\nu$. In the following part of the paper, we assume the absolute
continuity of $\mu,\nu$, and $\pi$ with respect to the Lebesgue measure for
well-defined entropy regularization.

##  3 Entropy-Regularized Optimal Transport between Multivariate Normal
Distributions

In this section, we provide a rigorous solution of entropy-regularized optimal
transport between two multivariate normal distributions. Throughout this
section, we adopt the squared Euclidean distance $\|x-y\|^{2}$ as the cost
function. To prove our theorem, we start by expressing $C_{\lambda}$ using
mean vectors and covariance matrices. The following lemma is a known result;
for example, see [18].

######  Lemma 3.1.

Let $X\sim P,Y\sim Q$ be two random variables on $\mathbb{R}^{n}$ with means
$\mu_{1},\mu_{2}$ and covariance matrices $\Sigma_{1},\Sigma_{2}$,
respectively. If $\pi(x,y)$ is a coupling of $P$ and $Q$, we have:

$\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}\|{x}-{y}\|^{2}\pi({x},{y})dxdy=\|\mu_{1}-\mu_{2}\|^{2}+\mathrm{tr}\left\\{\Sigma_{1}+\Sigma_{2}-2\mathrm{Cov}(X,Y)\right\\}.$
(12)

###### Proof.

Without loss of generality, we can assume $X$ and $Y$ are centralized,
because:

$\int\|(x-\mu_{1})-(y-\mu_{2})\|^{2}\pi(x,y)dxdy=\int\|x-y\|^{2}\pi(x,y)dxdy-\|\mu_{1}-\mu_{2}\|^{2}.$
(13)

Therefore, we have:

$\displaystyle\int\|{x}-{y}\|^{2}\pi({x},{y})d{x}d{y}$
$\displaystyle=E[\|X-Y\|^{2}]=E[\mathrm{tr}\\{(X-Y)(X-Y)^{\scalebox{0.6}{\rm
T}}\\}]$
$\displaystyle=\mathrm{tr}\left\\{\Sigma_{1}+\Sigma_{2}-2\mathrm{Cov}(X,Y)\right\\}.$
(14)

By adding $\|\mu_{1}-\mu_{2}\|^{2}$, we obtain (12). ∎

Lemma 3.1 shows that
$\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}\|{x}-{y}\|^{2}\pi({x},{y})d{x}d{y}$
can be parameterized by the covariance matrices
$\Sigma_{1},\Sigma_{2},\mathrm{Cov}(X,Y)$. Because $\Sigma_{1}$ and
$\Sigma_{2}$ are fixed, the infinite-dimensional optimization of the coupling
$\pi$ is a finite-dimensional optimization of covariance matrix
$\mathrm{Cov}(X,Y)$.

We prepare the following lemma to prove Theorem 3.3.

######  Lemma 3.2.

Under a fixed mean and covariance matrix, the probability measure that
maximizes the entropy is a multivariate normal distribution.

Lemma 3.2 is a particular case of the principle of maximum entropy [25], and
the proof can be found in [35] Theorem 3.1.

######  Theorem 3.3.

Let
$P\sim\mathcal{N}({\mu_{1}},\Sigma_{1}),Q\sim\mathcal{N}({\mu_{2}},\Sigma_{2})$
be two multivariate normal distributions. The optimal coupling $\pi$ of $P$
and $Q$ of the entropy-regularized optimal transport:

$C_{\lambda}(P,Q)=\inf_{\pi\in\Pi(P,Q)}\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}\|{x}-{y}\|^{2}\pi({x},{y})d{x}d{y}-4\lambda\mathrm{Ent}(\pi).\quad\quad\quad(\text{\textasteriskcentered})$
(15)

is expressed as:

$\pi\sim\mathcal{N}\left(\begin{pmatrix}{\mu}_{1}\\\ {\mu}_{2}\\\
\end{pmatrix},\begin{pmatrix}\Sigma_{1}&\Sigma_{\lambda}\\\
\Sigma_{\lambda}^{\scalebox{0.6}{\rm T}}&\Sigma_{2}\\\ \end{pmatrix}\right)$
(16)

where:

$\Sigma_{\lambda}:=\Sigma_{1}^{1/2}(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2}\Sigma_{1}^{-1/2}-{\lambda}I.$
(17)

Furthermore, $C_{\lambda}(P,Q)$ can be written as:

$\displaystyle C_{\lambda}(P,Q)=$
$\displaystyle\|\mu_{1}-\mu_{2}\|^{2}+\mathrm{tr}(\Sigma_{1}+\Sigma_{2}-2(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2})$
$\displaystyle-$ $\displaystyle
2{\lambda}\log|(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2}-\lambda
I|-2\lambda n\log(2\pi\lambda)-4\lambda n\log(2\pi)-2\lambda n$ (18)

and the relative entropy version can be written as:

$\tilde{C}_{\lambda}(P,Q)=C_{\lambda}(P,Q)+2{\lambda}\log|\Sigma_{1}||\Sigma_{2}|+4\lambda
n\\{\log(2\pi)+1\\}.$ (19)

We note that we use the regularization parameter $4\lambda$ in
$(\text{\textasteriskcentered})$ for the sake of simplicity.

###### Proof.

Although the first half of the proof can be derived directly from Lemma 3.2,
we provide a proof of this theorem by Lagrange calculus, which will be used
later for the extension to $q$-normal distributions. Now, we define an
optimization problem that is equivalent to the entropy-regularized optimal
transport as follows:

$\displaystyle\mathrm{minimize}\ $
$\displaystyle\int\|x-y\|^{2}\pi(x,y)dxdy-4\lambda\mathrm{Ent}(\pi)$ (20)
$\displaystyle\mathrm{subject\ to}\ \ $ $\displaystyle\int\pi(x,y)dx=q(y)\
\text{for}\ \forall\ y\in\mathbb{R}^{n},$ $\displaystyle\int\pi(x,y)dy=p(x)\
\text{for}\ \forall\ x\in\mathbb{R}^{n}\ .$ (21)

Here, $p(x)$ and $q(y)$ are probability density functions of $P$ and $Q$,
respectively. Let $\alpha(x)$, $\beta(y)$ be Lagrange multipliers that
correspond to the above two constraints. The Lagrangian function of (21) is
defined as:

$\displaystyle L(\pi,\alpha,\beta):$
$\displaystyle=\int\|x-y\|^{2}\pi(x,y)dxdy+4\lambda\int\pi(x,y)\log\pi(x,y)dxdy$
$\displaystyle-\int\alpha(x)\pi(x,y)dxdy+\int\alpha(x)p(x)dx$
$\displaystyle-\int\beta(y)\pi(x,y)dxdy+\int\beta(y)q(y)dy.$ (22)

Taking the functional derivative of (22) with respect to $\pi$, we obtain:

$\delta
L(\pi,\alpha,\beta)=\int\left(\|x-y\|^{2}+4\lambda\log\pi(x,y)-\alpha(x)-\beta(y)\right)\delta\pi(x,y)dxdy.$
(23)

By the fundamental lemma of the calculus of variations, we have:

$\pi({x},{y})\propto\exp\left(\alpha({x})+\beta({y})-\frac{\|{x}-{y}\|^{2}}{4\lambda}\right).$
(24)

Here, $\alpha(x),\beta(y)$ are determined from the constraints (21). We can
assume that $\pi$ is a $2n$-variate normal distribution, because for a fixed
covariance matrix $\mathrm{Cov}(X,Y)$, $-\mathrm{Ent}(\pi)$ takes the infimum
when the coupling $\pi$ is a multivariate normal distribution by Lemma 3.2.
Therefore, we can express $\pi$ by using $z=(x^{\scalebox{0.6}{\rm
T}},y^{\scalebox{0.6}{\rm T}})^{\scalebox{0.6}{\rm T}}$ and a covariance
matrix $\Sigma:=\mathrm{Cov}(X,Y)$ as:

$\pi(x,y)\propto\exp\left\\{-\frac{1}{2}z^{\scalebox{0.6}{\rm
T}}\begin{pmatrix}\Sigma_{1}&\Sigma_{\phantom{2}}\\\
\Sigma^{\scalebox{0.6}{\rm T}}&\Sigma_{2}\\\ \end{pmatrix}^{-1}z\right\\}.$
(25)

Putting:

$\begin{pmatrix}\tilde{\Sigma}_{1}&\tilde{\Sigma}\\\
\tilde{\Sigma}^{\scalebox{0.6}{\rm T}}&\tilde{\Sigma}_{2}\\\
\end{pmatrix}:=\begin{pmatrix}\Sigma_{1}&\Sigma_{\phantom{2}}\\\
\Sigma^{\scalebox{0.6}{\rm T}}&\Sigma_{2}\\\ \end{pmatrix}^{-1},$ (26)

we write:

$\displaystyle-\frac{1}{2}z^{\scalebox{0.6}{\rm
T}}\begin{pmatrix}\Sigma_{1}&\Sigma_{\phantom{2}}\\\
\Sigma^{\scalebox{0.6}{\rm T}}&\Sigma_{2}\\\ \end{pmatrix}^{-1}z$
$\displaystyle=-\frac{1}{2}\begin{pmatrix}x^{\scalebox{0.6}{\rm
T}}&y^{\scalebox{0.6}{\rm
T}}\end{pmatrix}\begin{pmatrix}\tilde{\Sigma}_{1}&\tilde{\Sigma}\\\
\tilde{\Sigma}^{\scalebox{0.6}{\rm T}}&\tilde{\Sigma}_{2}\\\
\end{pmatrix}\begin{pmatrix}x\\\ y\end{pmatrix}$ (27)
$\displaystyle=-\frac{1}{2}x^{\scalebox{0.6}{\rm
T}}\tilde{\Sigma}_{1}x-\frac{1}{2}y^{\scalebox{0.6}{\rm
T}}\tilde{\Sigma}_{2}y-x^{\scalebox{0.6}{\rm T}}\tilde{\Sigma}y.$ (28)

According to block matrix inversion formula [43],
$\tilde{\Sigma}=-\Sigma_{1}^{-1}\Sigma A^{-1}$ holds, where
$A:=\Sigma_{2}-\Sigma^{\scalebox{0.6}{\rm T}}\Sigma_{1}^{-1}\Sigma$ is
positive definite. Then, comparing the term $x^{\scalebox{0.6}{\rm T}}y$
between $\eqref{siki:lagmoto}$ and $\eqref{siki:tenkai}$, we obtain
$\Sigma_{1}^{-1}\Sigma A^{-1}=\frac{1}{2\lambda}I$ and:

$2\lambda\Sigma_{1}^{-1}\Sigma=A=\Sigma_{2}-\Sigma^{\scalebox{0.6}{\rm
T}}\Sigma_{1}^{-1}\Sigma.$ (29)

Here, $\Sigma_{1}^{-1}\Sigma=\Sigma^{\scalebox{0.6}{\rm T}}\Sigma_{1}^{-1}$
holds, because $A$ is a symmetric matrix, and thus, we obtain:

$\lambda\Sigma_{1}^{-1}\Sigma+\lambda\Sigma^{\scalebox{0.6}{\rm
T}}\Sigma_{1}^{-1}=\Sigma_{2}-\Sigma^{\scalebox{0.6}{\rm
T}}\Sigma_{1}^{-1}\Sigma.$ (30)

Completing the square of the above equation, we obtain:

$(\Sigma_{1}^{-1/2}(\Sigma+\lambda I)\Sigma_{1}^{1/2})^{\scalebox{0.6}{\rm
T}}(\Sigma_{1}^{-1/2}(\Sigma+\lambda
I)\Sigma_{1}^{1/2})=\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I$
(31)

Let $Q$ be an orthogonal matrix; then, (31) can be solved as:

$\Sigma_{1}^{-1/2}(\Sigma+\lambda
I)\Sigma_{1}^{1/2}=Q(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2}.$
(32)

We rearrange the above equation as follows:

$\Sigma_{1}^{1/2}(\Sigma_{1}^{-1}\Sigma)\Sigma_{1}^{1/2}+\lambda
I=Q(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2}.$ (33)

Because the left terms and
$(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2}$ are all
symmetric positive definite, we can conclude that $Q$ is the identity matrix
by the uniqueness of the polar decomposition. Finally, we obtain:

$\Sigma=\Sigma_{1}^{1/2}(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2}\Sigma_{1}^{-1/2}-{\lambda}I=:\Sigma_{\lambda}.$
(34)

We obtain (18) by the direct calculation of $C_{\lambda}$ using Lemma 3.1 with
this $\Sigma_{\lambda}$. ∎

The following corollary helps us to understand the properties of
$\Sigma_{\lambda}$.

######  Corollary 3.3.1.

Let $\nu_{\lambda,1}\leq\nu_{\lambda,2}\leq\dot{\leq}\nu_{\lambda,n}$ be the
eigenvalues of $\Sigma_{\lambda}$; then, $\nu_{\lambda,i}$ monotonically
decreases with $\lambda$ for any $i\in\\{1,2,\dot{,}n\\}$.

###### Proof.

Because
$\Sigma_{1}^{-1/2}\Sigma_{\lambda}\Sigma_{1}^{1/2}=(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2}-{\lambda}I$
has the same eigenvalues as $\Sigma_{\lambda}$, if we let $\\{\nu_{0,i}\\}$ be
the eigenvalues of $\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}$,
$\nu_{\lambda,i}=\sqrt{\nu_{0,i}+\lambda^{2}}-\lambda$, which is a
monotonically decreasing function of the regularization parameter $\lambda$. ∎

By the proof, for large $\lambda$, we can prove
$\Sigma_{1}^{-1/2}\Sigma_{\lambda}\Sigma_{1}^{1/2}\simeq\frac{1}{2\lambda}\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}$
by diagonalization and $\nu_{\lambda,i}\simeq\frac{1}{2\lambda}\nu_{0,i}$.
Thus, $\Sigma_{\lambda}\simeq\frac{1}{2\lambda}\Sigma_{1}\Sigma_{2}$, and each
element of $\Sigma_{\lambda}$ converges to zero as $\lambda\rightarrow\infty$.

We show how entropy regularization behaves in two simple experiments. We
calculate the entropy-regularized optimal transport cost
$\mathcal{N}\left(\tiny\begin{pmatrix}0\\\
0\end{pmatrix},\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}\normalsize\right)$ and
$\mathcal{N}\left(\tiny\begin{pmatrix}0\\\
0\end{pmatrix},\begin{pmatrix}2&-1\\\ -1&2\end{pmatrix}\normalsize\right)$ in
the original version and the relative entropy version in Figure 1. We separate
the entropy-regularized optimal transport cost into the transport cost term
and regularization term and display both of them.

Figure 1: Graph of the entropy-regularized optimal transport cost between
$\mathcal{N}\left(\tiny\begin{pmatrix}0\\\
0\end{pmatrix},\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}\normalsize\right)$ and
$\mathcal{N}\left(\tiny\begin{pmatrix}0\\\
0\end{pmatrix},\begin{pmatrix}2&-1\\\ -1&2\end{pmatrix}\normalsize\right)$
with respect to $\lambda$ from zero to $10$.

It is reasonable that as $\lambda\downarrow 0$, $\Sigma_{\lambda}$ converges
to
$\Sigma_{1}^{1/2}(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2})^{1/2}\Sigma_{1}^{-1/2}$,
which is equal to the original optimal coupling of nonregularized optimal
transport and as $\lambda\rightarrow\infty$, $\Sigma_{\lambda}$ converges to
$\boldsymbol{0}$. This is a special case of Corollary 3.3.1.The larger
$\lambda$ becomes, the less correlated the optimal coupling is. We visualize
this behavior by computing the optimal couplings of two one-dimensional normal
distributions in Figure 2.

Figure 2: Contours of the density functions of the entropy-regularized optimal
coupling of $\mathcal{N}(0,1)$ and $\mathcal{N}(5,2)$ in three different
parameters $\lambda=0.1,1,10$. All of the optimal couplings are two-variate
normal distributions.

The left panel shows the original version. The transport cost is always
positive, and the entropy regularization term can take both signs in general;
then, the sign and total cost depend on their balance. We note that the
transport cost as a function of $\lambda$ is bounded, whereas the entropy
regularization is not. The boundedness of the optimal cost is deduced from
(3.1) and Corollary 3.3.1, and the unboundedness of the entropy regularization
is due to the regularization parameter $\lambda$ multiplied by the entropy.
The right panel shows the relative entropy version. It always takes a non-
negative value. Furthermore, because the total cost is bounded by the value
for the independent joint distribution (which is always a feasible coupling),
both the transport cost and the relative entropy regularization regularization
term are also bounded. Nevertheless, the larger the regularization parameter
$\lambda$, the greater the influence of entropy regularization over the total
cost.

It is known that a specific Riemannian metric can be defined in the space of
multivariate normal distributions, which induces the Wasserstein distance
[51]. To understand the effect of entropy regularization, we illustrate how
entropy regularization deforms this geometric structure in Figure 3. Here, we
generate $100$ two-variate normal distributions
$\\{\mathcal{N}(0,\Sigma_{r,k})\\}_{r,k\in\\{1,2,\dot{,}10\\}}$, where
$\\{\Sigma_{r,k}\\}$ is defined as:

$\quad\Sigma_{r,k}=\begin{pmatrix}\cos\left(2\pi\cdot\frac{k}{10}\right)&-\sin\left(2\pi\cdot\frac{k}{10}\right)\\\
\sin\left(2\pi\cdot\frac{k}{10}\right)&\cos\left(2\pi\cdot\frac{k}{10}\right)\end{pmatrix}^{\scalebox{0.6}{\rm
T}}\begin{pmatrix}1&0\\\
0&\sqrt{\frac{r}{10}}\end{pmatrix}\begin{pmatrix}\cos\left(2\pi\cdot\frac{k}{10}\right)&-\sin\left(2\pi\cdot\frac{k}{10}\right)\\\
\sin\left(2\pi\cdot\frac{k}{10}\right)&\cos\left(2\pi\cdot\frac{k}{10}\right)\end{pmatrix}.$
(35)

To visualize the geometric structure of these two-variate normal
distributions, we compute the relative entropy-regularized optimal transport
cost $\tilde{C}_{\lambda}$ between each pairwise two-variate normal
distributions. Then, we apply multidimensional scaling [29] to embed them into
a plane (see Figure 3). We can see entropy regularization deforming the
geometric structure of the space of multivariate normal distributions. The
deformation for distributions close to the isotopic normal distribution is
more sensitive to the change in $\lambda$.

Figure 3: Multidimensional scaling of two-variate normal distributions. The
pairwise dissimilarities are given by the square root of the entropy-
regularized optimal transport cost $\tilde{C}_{\lambda}$ for three different
regularization parameters $\lambda=0,0.01,0.05$. Each ellipse in the figure
represents a contour of the density function
$\\{\mathcal{N}(0,\Sigma_{r,k})\\}$.

The following corollary states that if we allow orthogonal transformations of
two multivariate normal distributions with fixed covariance matrices, then the
minimum and maximum of $C_{\lambda}$ are attained when $\Sigma_{1}$ and
$\Sigma_{2}$ are diagonalizable by the same orthogonal matrix or,
equivalently, when the ellipsoidal contours of the two density functions are
aligned with the same orthogonal axes.

######  Corollary 3.3.2.

With the same settings as in Theorem 3.3, fix $\mu_{1}$, $\mu_{2}$,
$\Sigma_{1}$, and all eigenvalues of $\Sigma_{2}$. When $\Sigma_{1}$ is
diagonalized as $\Sigma_{1}=\Gamma^{\scalebox{0.6}{\rm
T}}\Lambda_{1}^{\downarrow}\Gamma$, where $\Lambda_{1}^{\downarrow}$ is the
diagonal matrix of the eigenvalues of $\Sigma_{1}$ in descending order and
$\Gamma$ is an orthogonal matrix,

  1. (i)

$C_{\lambda}(P,Q)$ is minimized by $\Sigma_{2}=\Gamma^{\scalebox{0.6}{\rm
T}}\Lambda_{2}^{\downarrow}\Gamma$ and

  2. (ii)

$C_{\lambda}(P,Q)$ is maximized by $\Sigma_{2}=\Gamma^{\scalebox{0.6}{\rm
T}}\Lambda_{2}^{\uparrow}\Gamma$,

where $\Lambda_{2}^{\downarrow}$ and $\Lambda_{2}^{\uparrow}$ are the diagonal
matrices of the eigenvalues of $\Sigma_{2}$ in descending and ascending order,
respectively. Therefore, neither the minimizer, nor the maximizer depend on
the choice of $\lambda$.

###### Proof.

Because $\mu_{1}$, $\mu_{2}$, $\Sigma_{1}$, and all eigenvalues of
$\Sigma_{2}$ are fixed,

$\displaystyle C_{\lambda}(P,Q)$ $\displaystyle=-2{\rm
tr}\left((\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2}\right)-\frac{\lambda}{2}\log|(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\lambda^{2}I)^{1/2}-\lambda
I|+({\rm constant})$ (36)
$\displaystyle=\sum_{i=1}^{n}-2(\nu_{i}+\lambda^{2})^{1/2}-\frac{\lambda}{2}\log\\{(\nu_{i}+\lambda^{2})^{1/2}-\lambda\\}+({\rm
constant})$ (37) $\displaystyle=\sum_{i=1}^{n}g_{\lambda}(\log(\nu_{i}))+({\rm
constant})$ (38)

where $\nu_{1}\leq\dots\leq\nu_{n}$ are the eigenvalues of
$\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}$ and:

$g_{\lambda}(x):=-2(e^{x}+\lambda^{2})^{1/2}-\frac{\lambda}{2}\log\\{(e^{x}+\lambda^{2})^{1/2}-\lambda\\}.$
(39)

Note that $g_{\lambda}(x)$ is a concave function, because:

$g_{\lambda}^{\prime\prime}(x)=-\frac{e^{x}(4e^{x}+7\lambda^{2})}{8(e^{x}+\lambda^{2})^{3/2}}<0.$
(40)

Let
$\nu_{1}^{\downarrow\downarrow}\leq\dots\leq\nu_{n}^{\downarrow\downarrow}$
and $\nu_{1}^{\downarrow\uparrow}\leq\dots\leq\nu_{n}^{\downarrow\uparrow}$ be
the eigenvalues of $\Lambda_{1}^{\downarrow}\Lambda_{2}^{\downarrow}$ and
$\Lambda_{1}^{\downarrow}\Lambda_{2}^{\uparrow}$, respectively. By Exercise
6.5.3 of [9] or Theorem 6.13 and Corollary 6.14 of [22],

$\displaystyle(\log(\nu_{i}^{\downarrow\uparrow}))\prec(\log(\nu_{i}))\prec(\log(\nu_{i}^{\downarrow\downarrow})),$
(41)

Here, for $(a_{i}),(b_{i})\in\mathbb{R}^{n}$ such that $a_{1}\geq\dots\geq
a_{n}$ and $b_{1}\geq\dots\geq b_{n}$, $(a_{i})\prec(b_{i})$ means:

$\sum_{i=1}^{k}a_{i}\leq\sum_{i=1}^{k}b_{i}\mbox{~{}~{}for~{}~{}}k=1,\dots,n-1,\mbox{~{}~{}and~{}~{}}\sum_{i=1}^{n}a_{i}=\sum_{i=1}^{n}b_{i}$
(42)

and $(a_{i})$ is said to be majorized by $(b_{i})$. Because $g_{\lambda}(x)$
is concave,

$g_{\lambda}(\log(\nu_{i}^{\downarrow\uparrow}))\prec^{\rm
w}g_{\lambda}(\log(\nu_{i}))\prec^{\rm
w}g_{\lambda}(\log(\nu_{i}^{\downarrow\downarrow})),$ (43)

where $\prec^{\rm w}$ represents weak supermajorization, i.e.,
$(a_{i})\prec^{\rm w}(b_{i})$ means:

$\sum_{i=k}^{n}a_{i}\geq\sum_{i=k}^{n}b_{i}\mbox{~{}~{}for~{}~{}}k=1,\dots,n$
(44)

(see Theorem 5.A.1 of [37], for example). Therefore,

$\sum_{i=1}^{n}g_{\lambda}(\log(\nu_{i}^{\downarrow\uparrow}))\geq\sum_{i=1}^{n}g_{\lambda}(\log(\nu_{i}))\geq\sum_{i=1}^{n}g_{\lambda}(\log(\nu_{i}^{\downarrow\downarrow})).$
(45)

As in Case (i) (or (ii)), the eigenvalues of
$\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}$ correspond to the eigenvalues of
$\Lambda_{1}^{\downarrow}\Lambda_{2}^{\downarrow}$ (or
$\Lambda_{1}^{\downarrow}\Lambda_{2}^{\uparrow}$, respectively), the corollary
follows. ∎

Note that a special case of Corollary 3.3.2 for the ordinary Wasserstein
metric ($\lambda=0$) has been studied in the context of fidelity and the Bures
distance in quantum information theory. See Lemma 3 of [36]. Their proof is
not directly applicable to our generalized result; thus, we used another
approach to prove it.

##  4 Extension to Tsallis Entropy Regularization

In this section, we consider a generalization of entropy-regularized optimal
transport. We now focus on the Tsallis entropy [52], which is a generalization
of the Shannon entropy and appears in nonequilibrium statistical mechanics. We
show that the optimal coupling of Tsallis entropy-regularized optimal
transport between two $q$-normal distributions is also a $q$-normal
distribution. We start by recalling the definition of the $q$-exponential
function and $q$-logarithmic function based on [52].

######  Definition 4.1.

Let $q$ be a real parameter, and let $u>0$. The $q$-logarithmic function is
defined as:

$\log_{q}(u):=\begin{cases}\frac{1}{1-q}(u^{1-q}-1)\quad\text{if}\ \ q\neq
1,\\\ \log(u)\quad\text{if}\ \ q=1\end{cases}$ (46)

and the $q$-exponential function is defined as:

$\exp_{q}(u):=\begin{cases}[1+(1-q)u]_{+}^{\frac{1}{1-q}}\quad\text{if}\ \
q\neq 1,\\\ \exp(u)\quad\text{if}\ \ q=1\end{cases}$ (47)

######  Definition 4.2.

Let $q<1$ or $1<q<1+\frac{2}{n}$; an $n$-variate $q$-normal distribution is
defined by two parameters: $\mu\in\mathbb{R}^{n}$ and a positive definite
matrix $\Sigma$, and its density function is:

$f(x):=\frac{1}{C_{q}(\Sigma)}\exp_{q}\left(-(x-\mu)^{\scalebox{0.6}{\rm
T}}\Sigma^{-1}(x-\mu)\right),$ (48)

where $C_{q}(\Sigma)$ is a normalizing constant. $\mu$ and $\Sigma$ are called
the location vector and scale matrix, respectively.

In the following, we write the multivariate $q$-normal distribution
$\mathcal{N}_{q}(\mu,\Sigma)$. We note that the property of the $q$-normal
distribution changes in accordance with $q$. The $q$-normal distribution has
an unbounded support for $1<q<\frac{2}{n}$ and a bounded support for $q<1$.
The second moment exists for $q<1+\frac{2}{n+2}$, and the covariance becomes
$\frac{1}{2+(n+2)(1-q)}\Sigma$. We remark that each $n$-variate
$\left(1+\frac{2}{\nu+n}\right)$-normal distribution is equivalent to an
$n$-variate $t$-distribution with $\nu$ degrees of freedom,

$\frac{\Gamma[(\nu+n)/2]}{\Gamma(\nu/2)\nu^{n/2}\pi^{n/2}|{\Sigma}|^{1/2}}\left[1+\frac{1}{\nu}({x}-{\mu})^{T}{\Sigma}^{-1}({x}-{\mu})\right]^{-(\nu+n)/2},$
(49)

for $1<q<1+\frac{2}{n+2}$ and an $n$-variate normal distribution for
$q\downarrow 1$.

######  Definition 4.3.

Let $p$ be a probability density function. The Tsallis entropy is defined as:

$S_{q}(p):=\int p(x)\log_{q}\frac{1}{p(x)}dx=\frac{1}{q-1}\left(1-\int
p(x)^{q}dx\right).$ (50)

Then, the Tsallis entropy-regularized optimal transport is defined as:

$\displaystyle\mathrm{minimize}\ $
$\displaystyle\int\|x-y\|^{2}\pi(x,y)dxdy-2\lambda S_{q}(\pi)$ (51)
$\displaystyle\mathrm{subject\ to}\ \ $ $\displaystyle\int\pi(x,y)dx=q(y)\
\text{for}\ \forall\ y\in\mathbb{R}^{n},$ $\displaystyle\int\pi(x,y)dy=p(x)\
\text{for}\ \forall\ x\in\mathbb{R}^{n}\ .$ (52)

The following lemma is a generalization of the maximum entropy principle for
the Shannon entropy shown in Section 2 of [14].

######  Lemma 4.4.

Let $P$ be a centered $n$-dimensional probability measure with a fixed
covariance matrix $\Sigma$; the maximizer of the Renyi $\alpha$-entropy:

$\frac{1}{1-\alpha}\log\int f(x)^{\alpha}\mathrm{d}x$ (53)

under the constraint is
$\mathcal{N}_{2-\alpha}(\mathbf{0},((n+2)\alpha-n)\Sigma)$ for
$\frac{n}{n+2}<\alpha<1$.

We note that the maximizers of the Renyi $\alpha$-entropy and the Tsallis
entropy with $q=\alpha$ coincide; thus, the above lemma also holds for the
Tsallis entropy. This is mentioned, for example, in Section 9 of [41].

To prove Theorem 4.6, we use the following property of multivariate
$t$-distributions, which is summarized in Chapter 1 of [28].

######  Lemma 4.5.

Let $X$ be a random vector following an $n$-variate $t$-distribution with
degree of freedom $\nu$. Considering a partition of the mean vector $\mu$ and
scale matrix $\Sigma$, such as:

${X}=\left(\begin{array}[]{l}{X}_{1}\\\
{X}_{2}\end{array}\right),\quad{\mu}=\left(\begin{array}[]{l}{\mu}_{1}\\\
{\mu}_{2}\end{array}\right),\quad\Sigma=\left(\begin{array}[]{ll}\Sigma_{11}&\Sigma_{12}\\\
\Sigma_{21}&\Sigma_{22}\end{array}\right),$ (54)

$X_{1}$ follows a $p$-variate $t$-distribution with degree of freedom $\nu$,
mean vector $\mu_{1}$, and scale matrix $\Sigma_{11}$, where $p$ is the
dimension of $X_{1}$.

Recalling the correspondence of the parameter of the multivariate $q$-normal
distribution and the degree of freedom of the multivariate $t$-distribution
$q=1+\frac{2}{\nu+n}$, we can obtain the following corollary.

######  Corollary 4.5.1.

Let $X$ be a random vector following an $n$-variate $q$-normal distribution
for $1<q<1+\frac{2}{n+2}$. Consider a partition of the mean vector $\mu$ and
scale matrix $\Sigma$ in the same way as in (54). Then, $X_{1}$ follows a
$p$-variate $\left(1+\frac{2(q-1)}{2-(n-p)(q-1)}\right)$-normal distribution
with mean vector $\mu_{1}$ and scale matrix $\Sigma_{11}$, where $p$ is the
dimension of $X_{1}$.

######  Theorem 4.6.

Let
$P\sim\mathcal{N}_{q}({\mu_{1}},\Sigma_{1}),Q\sim\mathcal{N}_{q}({\mu_{2}},\Sigma_{2})$
be $n$-variate $q$-normal distributions for $1<q<1+\frac{2}{n+2}$ and
$\tilde{q}=-\frac{2(q-1)}{2-n(q-1)}$; consider the Tsallis entropy-regularized
optimal transport:

$C_{\lambda}(P,Q)=\inf_{\pi\in\Pi(P,Q)}\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}\|{x}-{y}\|^{2}\pi({x},{y})d{x}d{y}-2\lambda
S_{1+\tilde{q}}(\pi).$ (55)

Then, there exists a unique
$\tilde{\lambda}=\tilde{\lambda}(q,\Sigma_{1},\Sigma_{2},\lambda)\in\mathbb{R}_{+}$
such that the optimal coupling $\pi$ of the entropy-regularized optimal
transport is expressed as:

$\pi\sim\mathcal{N}_{1-\tilde{q}}\left(\begin{pmatrix}{\mu}_{1}\\\
{\mu}_{2}\\\
\end{pmatrix},\begin{pmatrix}\Sigma_{1}&\Sigma_{\tilde{\lambda}}\\\
\Sigma_{\tilde{\lambda}}^{\scalebox{0.6}{\rm T}}&\Sigma_{2}\\\
\end{pmatrix}\right),$ (56)

where:

$\Sigma_{\tilde{\lambda}}:=\Sigma_{1}^{1/2}(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+{\tilde{\lambda}}^{2}I)^{1/2}\Sigma_{1}^{-1/2}-{\tilde{\lambda}}I.$
(57)

###### Proof.

The proof proceeds in a similar way as in Theorem 3.3. Let $\alpha\in L(P)$
and $\beta\in L(Q)$ be the Lagrangian multipliers. Then, the Lagrangian
function $L(\pi,\alpha,\beta)$ of (52) is defined as:

$\displaystyle L(\pi,\alpha,\beta):=$
$\displaystyle\int\|x-y\|^{2}\pi(x,y)dxdy-2\lambda\left\\{\frac{1}{\tilde{q}}\left(1-\int\pi(x,y)^{1+\tilde{q}}dxdy\right)\right\\}$
$\displaystyle-\int\alpha(x)\pi(x,y)dxdy+\int\alpha(x)p(x)dx$
$\displaystyle-\int\beta(y)\pi(x,y)dxdy+\int\beta(y)q(y)dy$ (58)

and the extremum of the Tsallis entropy-regularized optimal transport is
obtained by the functional derivative with respect to $\pi$,

$\pi(x,y)=\left(\frac{\tilde{q}}{{2(\tilde{q}+1)\lambda}}\left(-\alpha({x})-\beta({y})+\|{x}-{y}\|^{2}\right)\right)^{\frac{1}{\tilde{q}}}.$
(59)

Here, $\alpha$ and $\beta$ are quadratic polynomials by Lemma 4.4. To separate
the normalizing constant, we introduce a constant $c\in\mathbb{R}_{+}$, and
$\pi$ can be written as:

$\pi(x,y)=c^{\frac{1}{\tilde{q}}}\left(\tilde{\alpha}(x)+\tilde{\beta}({y})+\frac{\tilde{q}\|{x}-{y}\|^{2}}{2c(\tilde{q}+1)\lambda}\right)^{\frac{1}{\tilde{q}}},$
(60)

with quadratic functions $\tilde{\alpha}(x)$ and $\tilde{\beta}(y)$.

Let $\tilde{\lambda}=\frac{c(\tilde{q}+1)\lambda}{{\tilde{q}}}>0$. Then, by
the same argument as in the proof of Theorem 3.3 and using Corollary 4.5.1, we
obtain the scale matrix of $\pi$ as:

$\Sigma=\begin{pmatrix}\Sigma_{1}&\Sigma_{\tilde{\lambda}}\\\
\Sigma_{\tilde{\lambda}}^{\scalebox{0.6}{\rm T}}&\Sigma_{2}\\\ \end{pmatrix},$
(61)

where:

$\Sigma_{\tilde{\lambda}}=\Sigma_{1}^{1/2}(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}+\tilde{\lambda}^{2}I)^{1/2}\Sigma_{1}^{-1/2}-{\tilde{\lambda}}I.$
(62)

Let $z=(x^{\scalebox{0.6}{\rm T}},y^{\scalebox{0.6}{\rm
T}})^{\scalebox{0.6}{\rm T}}$ and $K_{\tilde{q}}=\int(1+z^{\scalebox{0.6}{\rm
T}}z)^{\frac{1}{\tilde{q}}}dz$; $\pi$ can be written as:

$\pi(x,y)=\frac{1}{K_{\tilde{q}}|\Sigma|}(1+z^{\scalebox{0.6}{\rm
T}}\Sigma^{-1}z)^{\frac{1}{\tilde{q}}}.$ (63)

The constant $c$ is determined by:

$\frac{1}{K_{\tilde{q}}|\Sigma|}=c^{\frac{1}{\tilde{q}}}.$ (64)

We will show that the above equation has a unique solution. Let
$\\{\tau\\}_{i=1}^{n}$ be the eigenvalues of
$(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2})^{1/2}$; $|\Sigma|$ can be
expressed as
$\prod_{i=1}^{2n}2\tilde{\lambda}(\sqrt{\tau_{i}^{2}+\tilde{\lambda}^{2}}-\tilde{\lambda})$.
We consider:

$\displaystyle f(c)$
$\displaystyle=\log(c^{\frac{1}{\tilde{q}}}{K_{\tilde{q}}|\Sigma|})$ (65)
$\displaystyle=\frac{1}{\tilde{q}}\log
c+\sum_{i=1}^{2n}\log(\sqrt{\tau_{i}^{2}+\tilde{\lambda}^{2}}-\tilde{\lambda})+2n\log(2\tilde{\lambda})+\log
K_{\tilde{q}}.$ (66)

Because $\tilde{q}<0$, $f(c)$ is a monotonic decreasing function, and
$\lim_{c\downarrow 0}f(c)=\infty$, $\lim_{c\to\infty}f(c)=-\infty$, (64) has a
unique positive solution, and $\tilde{\lambda}$ is determined uniquely. ∎

##  5 Entropy-Regularized Kantorovich Estimator

Many estimators are defined by minimizing the divergence or distance $\rho$
between probability measures, that is $\arg\min_{\mu}\rho(\mu,\nu)$ for a
fixed $\nu$. When $\rho$ is the Kullback–Leibler divergence, the estimator
corresponds to the maximum likelihood estimator. When $\rho$ is the
Wasserstein distance, the following estimator is called the minimum
Kantorovich estimator, according to [44]. In this section, we consider a
probability measure $Q^{*}$ that minimizes $C_{\lambda}(P,Q)$ for a fixed $P$
over $\mathcal{P}_{2}(\mathbb{R}^{n})$, the set of all probability measures on
$\mathbb{R}^{n}$ with finite second moment that are absolutely continuous with
respect to the Lebesgue measure. In other words, we define the entropy-
regularized Kantorovich estimator
$\arg\min_{Q\in\mathcal{P}_{2}(\mathbb{R}^{n})}C_{\lambda}(P,Q).$ The entropy-
regularized Kantorovich estimator for discrete probability measures was
studied in [5], Theorem 2. We obtain the entropy-regularized Kantorovich
estimator for continuous probability measures in the following theorem:

######  Theorem 5.1.

For a fixed $P\in\mathcal{P}_{2}(\mathbb{R}^{n})$,

$Q^{*}=\arg\min_{Q\in\mathcal{P}_{2}(\mathbb{R}^{n})}C_{\lambda}(P,Q)$ (67)

exists, and its density function can be written as:

$dQ^{*}=dP\star\phi_{\lambda},$ (68)

where $\phi_{\lambda}(x)$ is a density function of
$\mathcal{N}(0,\frac{\lambda}{2}I)$, and $\star$ denotes the convolution
operator.

We use the dual problem of the entropy-regularized optimal transport to prove
Theorem 5.1 (for details, see Proposition 2.1 of [7] or Section 3 of [12]).

######  Lemma 5.2.

The dual problem of entropy-regularized optimal transport can be written as:

$\displaystyle\mathcal{A}_{\lambda}(P,Q)=\sup_{\begin{subarray}{c}\alpha\in
L_{1}(P)\\\ \beta\in L_{1}(Q)\end{subarray}}$
$\displaystyle\int\alpha(x)p(x)dx+\int\beta(y)q(y)dy$
$\displaystyle-\lambda\int\exp\left\\{\frac{\alpha(x)+\beta(y)-\|x-y\|^{2}}{\lambda}\right\\}dxdy.$
(69)

Moreover, $\mathcal{A}_{\lambda}(P,Q)=C_{\lambda}(P,Q)$ holds.

Now, we prove Theorem 5.1.

###### Proof.

Let $Q^{*}$ be the minimizer of $\min_{Q}C_{\lambda}(P,Q)$. Applying Lemma
5.2, there exist $\alpha^{*}\in L_{1}(P)$ and $\beta^{*}\in L_{1}(Q^{*})$ such
that:

$\displaystyle C_{\lambda}(P,Q^{*})=\mathcal{A}_{\lambda}(P,Q^{*})$
$\displaystyle=\int\alpha^{*}(x)p(x)dx+\int\beta^{*}(y)q^{*}(y)dy$
$\displaystyle-\lambda\int\exp\left\\{\frac{\alpha^{*}(x)+\beta^{*}(y)-\|x-y\|^{2}}{\lambda}\right\\}dxdy.$
(70)

Now, $\mathcal{A}_{\lambda}(P,Q^{*})$ is the minimum value of
$\mathcal{A}_{\lambda}$, such that the variation
$\delta\mathcal{A}_{\lambda}(P,Q^{*})$ is always zero. Then,

$\delta\mathcal{A}_{\lambda}(P,Q^{*})=\int\beta^{*}(y)\delta
q^{*}(y)dy=0\Rightarrow\beta^{*}\equiv 0$ (71)

holds, and the optimal coupling of $P,Q$ can be written as:

$\displaystyle\pi^{*}(x,y)$
$\displaystyle=\exp\left\\{\frac{\alpha^{*}(x)+\beta^{*}(y)}{\lambda}-\frac{\|x-y\|^{2}}{\lambda}\right\\}$
(72)
$\displaystyle=\exp\left\\{\frac{\alpha^{*}(x)}{\lambda}\right\\}\exp\left\\{-\frac{\|x-y\|^{2}}{\lambda}\right\\}.$
(73)

Moreover, we can obtain a closed-form of $\alpha^{*}(x)$ as follows from the
equation $\int\pi(x,y)dy$ $=p(x)$:

$\frac{\alpha^{*}(x)}{\lambda}=\log
p(x)-\log\int\exp\left\\{-\frac{\|x-y\|^{2}}{\lambda}\right\\}dy=\log
p(x)-\frac{n}{2}\log(\pi\lambda).$ (74)

Then, by calculating the marginal distribution of $\pi(x,y)$ with respect to
$x$, we can obtain:

$q^{*}(y)=\int\frac{1}{(\pi{\lambda})^{\frac{n}{2}}}\exp\left\\{-\frac{\|x-y\|^{2}}{\lambda}\right\\}p(x)dx=(p\star\phi_{\lambda})(y).$
(75)

Therefore, we conclude that a probability measure $Q$ that minimizes
$C_{\lambda}(P,Q)$ is expressed as (75). ∎

It should be noted that when $P$ in Theorem 5.1 are multivariate normal
distributions, $Q^{*}$ and $P$ are simultaneously diagonalizable by a direct
consequence of the theorem. This is consistent with the result of Corollary
3.3.2(1) for minimization when all eigenvalues are fixed.

We can determine that the entropy-regularized Kantorovich estimator is a
measure convolved with an isotropic multivariate normal distribution scaled by
the regularization parameter $\lambda$. This is similar to the idea of prior
distributions in the context of Bayesian inference. Applying Theorem 5.1, the
entropy-regularized Kantorovich estimator of the multivariate normal
distribution $\mathcal{N}(\mu,\Sigma)$ is
$\mathcal{N}(\mu,\Sigma+\frac{\lambda}{2}I)$.

##  6 Numerical Experiments

In this section, we introduce experiments that show the statistical efficiency
of entropy regularization in Gaussian settings. We consider two different
setups, estimating covariance matrices (Section 6.1) and the entropy-
regularized Wasserstein barycenter (Section 6.2). To obtain the entropy-
regularized Wasserstein barycenter, we adopt the Newton–Schulz method and a
manifold optimization method, which are explained in Sections 6.3 and 6.4,
respectively.

###  6.1 Estimation of Covariance Matrices

We provide a covariance estimation method based on entropy-regularized optimal
transport. Let $P=\mathcal{N}(\mu,\Sigma)$ be an $n$-variate normal
distribution. We define an entropy-regularized Kantorovich estimator
$\hat{P}_{\lambda}$, that is,

$\hat{P}_{\lambda}=\arg\min_{Q}C_{\lambda}(P,Q).$ (76)

We generate some samples from $\mathcal{N}(\mu,\Sigma)$ and estimate the mean
and covariance matrix. We compare the maximum likelihood estimator
$\hat{P}_{\mathrm{MLE}}=\mathcal{N}(\hat{\mu}_{\mathrm{MLE}},\hat{\Sigma}_{\mathrm{MLE}})$
and $\hat{P}_{\lambda}$ with respect to the prediction error:

$\mathrm{KL}(P,\hat{P}_{\mathrm{MLE}}),\ \ \mathrm{KL}(P,\hat{P}_{\lambda}).$
(77)

In our experiment, the dimension $n$ is set to $5,15,30$, and the sample size
is set to $60,120$. The experiment proceeds as follows.

  1. 1.

Obtain a random sample of size 60 (or 120) from $\mathcal{N}(0,\Sigma)$ and
its sample covariance matrix $\hat{\Sigma}$.

  2. 2.

Obtain the entropy-regularized minimum Kantorovich estimator of $\hat{\Sigma}$
obtained in Step 1.

  3. 3.

Compute the prediction error between $\Sigma$ and the entropy-regularized
minimum Kantorovich estimator of $\hat{\Sigma}$

  4. 4.

Repeat the above steps 1000 times and obtain a confidence interval of the
prediction error.

Table 1 shows the average prediction error of the MLE and entropy-regularized
Kantorovich estimator of covariance matrices from 60 samples from an n-variate
normal distribution with the 95% confidential interval. We can see that the
prediction error is smaller than the maximum likelihood estimator under
adequately small $\lambda$ for $n=15,30$, but not for $n=5$. Moreover, the
decrease in the prediction error is larger for $n=30$ than for $n=15$, which
indicates that the entropy regularization is effective in a high dimension. On
the other hand, Table 2 shows in all cases that the decreases in the
prediction error are more moderate than Table 1. We can see that this is due
to the increase in the sample size. Then, we can conclude that the entropy
regularization is effective in a high-dimensional setting with a small sample
size.

Table 1: Average prediction error of the MLE and entropy-regularized Kantorovich estimator of covariance matrices from 60 samples from an $n$-variate normal distribution with the 95% confidential interval. $\lambda$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=5$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=15$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=30$
---|---|---|---
0(MLE) | 0.062 $\pm\ 0.005$ | 1.346 $\pm\ 0.022$ | 10.69 $\pm\ 0.112$
0.01 | 0.051 $\pm\ 0.005$ | 1.242 $\pm\ 0.021$ | 8.973 $\pm\ 0.087$
0.1 | 0.104 $\pm\ 0.004$ | 0.841 $\pm\ 0.013$ | 4.180 $\pm\ 0.033$
0.5 | 0.647 $\pm\ 0.003$ | 0.931 $\pm\ 0.007$ | 3.093 $\pm\ 0.010$
1.0 | 1.166 $\pm\ 0.003$ | 1.670 $\pm\ 0.006$ | 5.075 $\pm\ 0.009$
Table 2: Average prediction error of the MLE and entropy-regularized Kantorovich estimator of covariance matrices from 120 samples from an $n$-variate normal distribution with the 95% confidential interval. $\lambda$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=5$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=15$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=30$
---|---|---|---
0(MLE) | 0.024 $\pm\ 0.002$ | 0.490 $\pm\ 0.007$ | 2.810 $\pm\ 0.021$
0.01 | 0.020 $\pm\ 0.002$ | 0.459 $\pm\ 0.006$ | 2.528 $\pm\ 0.018$
0.1 | 0.101 $\pm\ 0.002$ | 0.397 $\pm\ 0.005$ | 1.700 $\pm\ 0.001$
0.5 | 0.659 $\pm\ 0.002$ | 0.875 $\pm\ 0.004$ | 2.833 $\pm\ 0.005$
1.0 | 1.180 $\pm\ 0.002$ | 1.730 $\pm\ 0.004$ | 5.124 $\pm\ 0.005$

###  6.2 Estimation of the Wasserstein Barycenter

A barycenter with respect to the Wasserstein distance is definable [1] and is
widely used for image interpolation and 3D object interpolation tasks with
entropy regularization [5, 49].

######  Definition 6.1.

Let $\\{Q_{i}\\}_{i=1}^{m}$ be a set of probability measures in
$\mathcal{P}(\mathbb{R}^{n})$. The barycenter with respect to $C_{\lambda}$
(entropy-regularized Wasserstein barycenter) is defined as:

$\arg\min_{P\in\mathcal{P}(\mathbb{R}^{n})}\sum_{i=1}^{m}C_{\lambda}(P,Q_{i}).$
(78)

Now, we restrict $P$ and $\\{Q_{i}\\}_{i=1}^{m}$ to be multivariate normal
distributions and apply our theorem to illustrate the effect of entropy
regularization.

The experiment proceeds as follows. The dimensionality and the sample size
were set the same as in the experiments in Section 6.1.

  1. 1.

Obtain a random sample of size 60 (or 120) from $\mathcal{N}(0,\Sigma)$ and
its sample covariance matrix $\hat{\Sigma}$.

  2. 2.

Repeat Step 1 three times, and obtain $\\{\hat{\Sigma}\\}_{i=1}^{3}$.

  3. 3.

Obtain the barycenter of $\\{\hat{\Sigma}_{i}\\}_{i=1}^{3}$.

  4. 4.

Compute the prediction error between $\Sigma$ and the barycenter obtained in
step 3.

  5. 5.

Repeat the above steps 100 times and obtain a confidence interval of the
prediction error.

We show the results for several values of the regularization parameter
$\lambda$ in Tables 3 and 4. A decrease in the prediction error can be seen in
Table 3 for $n=30$, as well as Tables 1 and 2. However, because the
computation of the entropy-regularized Wasserstein barycenter uses more data
than that of the minimum Kantorovich estimator, the decrease in the prediction
error is mild. The entropy-regularized Kantorovich estimator is a special case
of the entropy-regularized Wasserstein barycenter (78) for $m=1$. Our
experiments show that the appropriate range of $\lambda$ to decrease the
prediction error depends on $m$ and becomes narrow as $m$ increases. In
addition, we note that there is a small decrease in the prediction error in
Table 4 for $n=30$.

Table 3: Average prediction error of the entropy-regularized barycenter with the 95% confidential interval (random sample of size 60). $\lambda$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=5$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=15$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=30$
---|---|---|---
0 | 0.455 $\pm\ 0.395$ | 1.318 $\pm\ 0.006$ | 4.875 $\pm\ 0.035$
0.001 | 0.429 $\pm\ 0.396$ | 1.318 $\pm\ 0.004$ | 4.887 $\pm\ 0.036$
0.01 | 0.434 $\pm\ 0.270$ | 1.344 $\pm\ 0.006$ | 4.551 $\pm\ 0.164$
0.025 | 0.780 $\pm\ 0.223$ | 1.456 $\pm\ 0.064$ | 5.710 $\pm\ 0.536$
0.005 | 1.047 $\pm\ 0.029$ | 1.537 $\pm\ 0.064$ | 7.570 $\pm\ 0.772$
Table 4: Average prediction error of the entropy-regularized barycenter with the 95% confidential interval (random sample of size 120). $\lambda$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=5$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=15$ | $\mathrm{KL}(P,\hat{P}_{\mathrm{W}}),n=30$
---|---|---|---
0 | 0.154 $\pm\ 0.600$ | 1.303 $\pm\ 0.010$ | 5.091 $\pm\ 0.035$
0.001 | 0.212 $\pm\ 0.070$ | 1.305 $\pm\ 0.010$ | 5.072 $\pm\ 0.037$
0.01 | 0.306 $\pm\ 0.046$ | 1.328 $\pm\ 0.008$ | 5.274 $\pm\ 0.252$
0.025 | 0.671 $\pm\ 0.028$ | 1.337 $\pm\ 0.073$ | 5.851 $\pm\ 0.424$
0.005 | 1.109 $\pm\ 0.063$ | 1.603 $\pm\ 0.184$ | 8.072 $\pm\ 0.725$

###  6.3 Gradient Descent on $\mathrm{Sym}_{+}(n)$

We use a gradient descent method to compute the entropy-regularized
barycenter. Applying the gradient descent method to the loss function defined
by the Wasserstein distance was proposed in [40]. This idea is extendable to
entropy-regularized optimal transport. The detailed algorithm is shown below.
Because $C_{\lambda}(P,Q)$ is a function of a positive definite matrix, we
used a manifold gradient descent algorithm on the manifold of positive
definite matrices.

We review the manifold gradient descent algorithm used in our numerical
experiment. Let $\mathrm{Sym}_{+}(n)$ be the manifold of $n$-dimensional
positive definite matrices. We require a formula for a gradient operator and
the inner product of $\mathrm{Sym}_{+}(n)$ in the gradient descent algorithm.
In this paper, we use the following inner product from [9], Chapter 6. For a
fixed $X\in\mathrm{int}(\mathrm{Sym}_{+}(n))$, we define an inner product of
$\mathrm{Sym}_{+}(n)$ as:

$g_{X}(Y,Z)=\operatorname{tr}\left(YX^{-1}ZX^{-1}\right),\
Y,Z\in\mathrm{Sym}_{+}(n),$ (79)

Equation (79) is the best choice in terms of the convergence speed according
to [26]. Let $f:\mathrm{Sym}_{+}(n)\rightarrow\mathbb{R}$ be a differential
matrix function. Then, the induced gradient of $f$ under (79) is:

$\operatorname{grad}f(X)=X\left(\frac{\partial f(X)}{\partial X}\right)X.$
(80)

We consider the updating step after obtaining the gradient of $f$.
$\mathrm{grad}f(X)$ is an element of the tangent space, and we have to project
it to $\mathrm{Sym}_{+}(n)$. This projection map is called a retraction. It is
known that the Riemannian metric $g_{X}$ leads to the following retraction:

$\exp_{X}x=X\operatorname{Exp}\left(X^{-1}x\right)$, where $\mathrm{Exp}$ is
the matrix exponential. Then, the corresponding gradient descent method
becomes as shown in Algorithm 1.

###  6.4 Approximate the Matrix Square Root

To compute the gradient of the square root of a matrix in the objective
function, we approximate it using the Newton–Schulz method [23], which can be
implemented by matrix operations as shown in Algorithm 2. It is amenable to
automatic differentiation, such that we can easily apply the gradient descent
method to our algorithm.

Algorithm 1 Gradient descent on the manifold of positive definite matrices.

0: $f(X)$

initialize $X$

while no convergence do

$\eta:\ $step size

$\mathrm{grad}\leftarrow X\left(\frac{\partial f(X)}{\partial X}\right)X$

$X\leftarrow\exp_{X}(\eta\mathrm{grad})=X\mathrm{Exp}(\eta
X^{-1}\mathrm{grad})$

end while

$X$

Algorithm 2 Newton–Schulz method.

0: $A\in\mathrm{Sym}_{+}(n),\epsilon>0$

$Y\leftarrow\frac{A}{(1+\epsilon)\|A\|},\ \ Z\leftarrow I$

while no convergence do

$T\leftarrow(3I-ZY)/2$

$Y\leftarrow YT,\ Z\leftarrow TZ$

end while

$\sqrt{(1+\epsilon)\|A\|}Y$

##  7 Conclusions and Future Work

In this paper, we studied entropy-regularized optimal transport and derived
several result. We summarize these as follows and add notes on future work.

  * •

We obtain the explicit form of entropy-regularized optimal transport between
two multivariate normal distributions and derived Corollaries 3.3.1 and 3.3.2,
which clarified the properties of optimal coupling. Furthermore, we
demonstrate experimentally how entropy regularization affects the Wasserstein
distance, the optimal coupling, and the geometric structure of multivariate
normal distributions. Overall, the properties of optimal coupling were
revealed both theoretically and experimentally. We expect that the explicit
formula can be a replacement for the existing methodology using the
(nonregularized) Wasserstein distance between normal distributions (for
example, [40, 21]).

  * •

Theorem 4.6 derives the explicit form of the optimal coupling of the Tsallis
entropy-regularized optimal transport between multivariate $q$-normal
distributions. The optimal coupling of the Tsallis entropy-regularized optimal
transport between multivariate $q$-normal distributions is also a multivariate
$q$-normal distribution, and the obtained result has an analogy to that of the
normal distribution. We believe that this result can be extended to other
elliptical distribution families.

  * •

The entropy-regularized Kantorovich estimator of a probability measure in
$\mathcal{P}_{2}(\mathbb{R})$ is the convolution of a multivariate normal
distribution and its own density function. Our experiments show that both the
entropy-regularized Kantorovich estimator and the Wasserstein barycenter of
multivariate normal distributions outperform the maximum likelihood estimator
in the prediction error for adequately selected $\lambda$ in a high
dimensionality and small sample setting. As future work, we want to show the
efficiency of entropy regularization using real data.

##  8 acknowledgements

This work was supported by RIKEN AIP and JSPS KAKENHI (JP19K03642,
JP19K00912).

## References

  * [1] Martial Agueh and Guillaume Carlier.  Barycenters in the Wasserstein space.  SIAM Journal on Mathematical Analysis, 43(2):904–924, 2011.
  * [2] Zeyuan Allen-Zhu, Yuanzhi Li, Rafael Oliveira, and Avi Wigderson.  Much faster algorithms for matrix scaling.  In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 890–901. IEEE, 2017.
  * [3] Jason Altschuler, Jonathan Weed, and Philippe Rigollet.  Near-linear time approximation algorithms for optimal transport via sinkhorn iteration.  arXiv preprint arXiv:1705.09634, 2017.
  * [4] Shun-ichi Amari, Ryo Karakida, and Masafumi Oizumi.  Information geometry connecting Wasserstein distance and Kullback–Leibler divergence via the entropy-relaxed transportation problem.  Information Geometry, 1(1):13–37, 2018.
  * [5] Shun-ichi Amari, Ryo Karakida, Masafumi Oizumi, and Marco Cuturi.  Information geometry for regularized optimal transport and barycenters of patterns.  Neural computation, 31(5):827–848, 2019.
  * [6] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 214–223. PMLR, 2017.
  * [7] Genevay Aude, Marco Cuturi, Gabriel Peyré, and Francis Bach.  Stochastic optimization for large-scale optimal transport.  arXiv preprint arXiv:1605.08527, 2016.
  * [8] Yogesh Balaji, Hamed Hassani, Rama Chellappa, and Soheil Feizi.  Entropic GANs meet VAEs: A statistical approach to compute sample likelihoods in GANs.  In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 414–423. PMLR, 09–15 Jun 2019.
  * [9] Rajendra Bhatia.  Positive definite matrices, volume 24.  Princeton university press, 2009.
  * [10] Jose Blanchet, Arun Jambulapati, Carson Kent, and Aaron Sidford.  Towards optimal running times for optimal transport.  arXiv preprint arXiv:1810.07717, 2018.
  * [11] Mathieu Blondel, Vivien Seguy, and Antoine Rolet.  Smooth and sparse optimal transport.  In International Conference on Artificial Intelligence and Statistics, pages 880–889. PMLR, 2018.
  * [12] Christian Clason, Dirk A Lorenz, Hinrich Mahler, and Benedikt Wirth.  Entropic regularization of continuous optimal transport problems.  Journal of Mathematical Analysis and Applications, 494(1):124432, 2021.
  * [13] Michael B Cohen, Aleksander Madry, Dimitris Tsipras, and Adrian Vladu.  Matrix scaling and balancing via box constrained newton’s method and interior point methods.  In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 902–913. IEEE, 2017.
  * [14] Jose Costa, Alfred Hero, and Christophe Vignat.  On solutions to multivariate maximum $\alpha$-entropy problems.  In International Workshop on Energy Minimization Methods in Computer Vision and Pattern Recognition, pages 211–226. Springer, 2003.
  * [15] Nicolas Courty, Rémi Flamary, Devis Tuia, and Alain Rakotomamonjy.  Optimal transport for domain adaptation.  IEEE transactions on pattern analysis and machine intelligence, 39(9):1853–1865, 2016.
  * [16] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In Advances in neural information processing systems, pages 2292–2300, 2013.
  * [17] Marco Cuturi and Gabriel Peyré.  A smoothed dual approach for variational wasserstein problems.  SIAM Journal on Imaging Sciences, 9(1):320–343, 2016.
  * [18] DC Dowson and BV Landau.  The Fréchet distance between multivariate normal distributions.  Journal of multivariate analysis, 12(3):450–455, 1982.
  * [19] Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin.  Computational optimal transport: Complexity by accelerated gradient descent is better than by sinkhorn’s algorithm.  In International conference on machine learning, pages 1367–1376. PMLR, 2018.
  * [20] Charlie Frogner, Chiyuan Zhang, Hossein Mobahi, Mauricio Araya, and Tomaso A Poggio.  Learning with a Wasserstein loss.  In Advances in Neural Information Processing Systems, pages 2053–2061, 2015.
  * [21] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [22] Fumio Hiai and Dénes Petz.  Introduction to matrix analysis and applications.  Springer Science & Business Media, 2014.
  * [23] Nicholas J Higham.  Newton’s method for the matrix square root.  Mathematics of Computation, 46(174):537–549, 1986.
  * [24] Hicham Janati, Boris Muzellec, Gabriel Peyré, and Marco Cuturi.  Entropic optimal transport between (unbalanced) gaussian measures has a closed form.  In Advances in neural information processing systems, 2020.
  * [25] Edwin T Jaynes.  Information theory and statistical mechanics.  Physical review, 106(4):620, 1957.
  * [26] Ben Jeuris, Raf Vandebril, and Bart Vandereycken.  A survey and comparison of contemporary algorithms for computing the matrix geometric mean.  Electronic Transactions on Numerical Analysis, 39(ARTICLE):379–402, 2012.
  * [27] Leonid V Kantorovich.  On the translocation of masses.  In Dokl. Akad. Nauk. USSR (NS), volume 37, pages 199–201, 1942\.
  * [28] Samuel Kotz and Saralees Nadarajah.  Multivariate t-distributions and their applications.  Cambridge University Press, 2004.
  * [29] Joseph B Kruskal.  Nonmetric multidimensional scaling: a numerical method.  Psychometrika, 29(2):115–129, 1964.
  * [30] Yin Tat Lee and Aaron Sidford.  Path finding methods for linear programming: Solving linear programs in o (vrank) iterations and faster algorithms for maximum flow.  In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pages 424–433. IEEE, 2014.
  * [31] Jing Lei et al.  Convergence and concentration of empirical measures under Wasserstein distance in unbounded functional spaces.  Bernoulli, 26(1):767–798, 2020.
  * [32] Tianyi Lin, Nhat Ho, and Michael Jordan.  On efficient optimal transport: An analysis of greedy and accelerated mirror descent algorithms.  In International Conference on Machine Learning, pages 3982–3991. PMLR, 2019.
  * [33] Tianyi Lin, Nhat Ho, and Michael I Jordan.  On the efficiency of the sinkhorn and greenkhorn algorithms and their acceleration for optimal transport.  arXiv preprint arXiv:1906.01437, 2019.
  * [34] Anton Mallasto, Augusto Gerolin, and Hà Quang Minh.  Entropy-regularized $2$-wasserstein distance between gaussian measures.  arXiv preprint arXiv:2006.03416, 2020.
  * [35] Kanti V Mardia.  Characterizations of directional distributions.  In A Modern Course on Statistical Distributions in Scientific Work, pages 365–385. Springer, 1975.
  * [36] Damian Markham, Jarosław Adam Miszczak, Zbigniew Puchała, and Karol Życzkowski.  Quantum state discrimination: a geometric approach.  Physical Review A, 77(4):042–111, 2008.
  * [37] Albert W Marshall, Ingram Olkin, and Barry C Arnold.  Inequalities: theory of majorization and its applications, volume 143.  Springer, 1979.
  * [38] Gonzalo Mena and Jonathan Niles-Weed.  Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem.  In Advances in Neural Information Processing Systems, pages 4543–4553, 2019.
  * [39] Gaspard Monge.  Mémoire sur la théorie des déblais et des remblais.  Histoire de l’Académie Royale des Sciences de Paris, 1781.
  * [40] Boris Muzellec and Marco Cuturi.  Generalizing point embeddings using the Wasserstein space of elliptical distributions.  In Advances in Neural Information Processing Systems, pages 10237–10248, 2018.
  * [41] Jan Naudts.  Generalised thermostatistics.  Springer Science & Business Media, 2011.
  * [42] Atsushi Nitanda and Taiji Suzuki.  Gradient layer: Enhancing the convergence of adversarial training for generative models.  In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 1008–1016. PMLR, 2018.
  * [43] KB Petersen, MS Pedersen, et al.  The matrix Cookbook.  Technical University of Denmark, 15, 2008.
  * [44] Gabriel Peyré, Marco Cuturi, et al.  Computational optimal transport.  Foundations and Trends® in Machine Learning, 11(5-6):355–607, 2019.
  * [45] Kent Quanrud.  Approximating optimal transport with linear programs.  arXiv preprint arXiv:1810.05957, 2018.
  * [46] Philippe Rigollet and Jonathan Weed.  Entropic optimal transport is maximum-likelihood deconvolution.  Comptes Rendus Mathematique, 356(11-12):1228–1235, 2018.
  * [47] Yossi Rubner, Carlo Tomasi, and Leonidas J Guibas.  The Earth Mover’s Distance as a metric for Image Retrieval.  International journal of computer vision, 40(2):99–121, 2000.
  * [48] Richard Sinkhorn and Paul Knopp.  Concerning nonnegative matrices and doubly stochastic matrices.  Pacific Journal of Mathematics, 21(2):343–348, 1967.
  * [49] Justin Solomon, Fernando De Goes, Gabriel Peyré, Marco Cuturi, Adrian Butscher, Andy Nguyen, Tao Du, and Leonidas Guibas.  Convolutional Wasserstein Distances: Efficient optimal transportation on geometric domains.  ACM Transactions on Graphics (TOG), 34(4):1–11, 2015.
  * [50] Sho Sonoda and Noboru Murata.  Transportation analysis of denoising autoencoders: A novel method for analyzing deep neural networks.  arXiv preprint arXiv:1712.04145, 2017.
  * [51] Asuka Takatsu et al.  Wasserstein geometry of Gaussian measures.  Osaka Journal of Mathematics, 48(4):1005–1026, 2011.
  * [52] Constantino Tsallis.  Possible generalization of Boltzmann-Gibbs statistics.  Journal of statistical physics, 52(1-2):479–487, 1988.
  * [53] Cédric Villani.  Optimal transport: old and new, volume 338.  Springer Science & Business Media, 2008.
  * [54] Pascal Vincent, Hugo Larochelle, Yoshua Bengio, and Pierre-Antoine Manzagol.  Extracting and composing robust features with denoising autoencoders.  In Proceedings of the 25th international conference on Machine learning, pages 1096–1103, 2008.