# Unsupervised Ground Metric Learning
using Wasserstein Eigenvectors

Geert-Jan Huizing†111Computational Systems Biology Team, Institut de Biologie
de l’Ecole Normale Supérieure, CNRS, INSERM, Ecole Normale Supérieure,
Université PSL, 75005, Paris, France , Laura Cantini∗, Gabriel
Peyré222Département de mathématiques et applications de l’Ecole Normale
Supérieure, CNRS, Ecole Normale Supérieure, Université PSL, 75005, Paris,
France CNRS and ENS, PSL University<EMAIL_ADDRESS>
###### Abstract

Optimal Transport (OT) defines geometrically meaningful “Wasserstein”
distances, used in machine learning applications to compare probability
distributions. However, a key bottleneck is the design of a “ground” cost
which should be adapted to the task under study. In most cases, supervised
metric learning is not accessible, and one usually resorts to some ad-hoc
approach. Unsupervised metric learning is thus a fundamental problem to enable
data-driven applications of Optimal Transport. In this paper, we propose for
the first time a canonical answer by computing the ground cost as a positive
eigenvector of the function mapping a cost to the pairwise OT distances
between the inputs. This map is homogeneous and monotone, thus framing
unsupervised metric learning as a non-linear Perron-Frobenius problem. We
provide criteria to ensure the existence and uniqueness of this eigenvector.
In addition, we introduce a scalable computational method using entropic
regularization, which – in the large regularization limit – operates a
principal component analysis dimensionality reduction. We showcase this method
on synthetic examples and datasets. Finally, we apply it in the context of
biology to the analysis of a high-throughput single-cell RNA sequencing
(scRNA-seq) dataset, to improve cell clustering and infer the relationships
between genes in an unsupervised way.

##  1 Introduction

Choosing an appropriate metric is a core ingredient for visualization and
learning on a high dimensional data set. This choice can be guided or improved
by leveraging metric learning methods [5]. This prior step is even more
important when making use of Optimal Transport (OT) methods. OT indeed lifts a
“ground” metric into a distance between discrete histograms or more general
probability distributions.

In this work we propose for the first time an unsupervised ground metric
learning method for OT. While existing approaches operate in a supervised way,
our method is similar in spirit to principal component analysis, which
formulates the problem as an eigenvector (or more generally singular vector)
computation problem.

###  1.1 Previous Works

#### Optimal Transport

While the initial proposal of Monge [42] formulates the OT problem as an
optimal matching problem, its modern and tractable formulation by [30] is a
linear program detailed in Section 3.1. Besides its use to define matchings
and couplings between distributions, the main feature of OT is that the
transportation value induces a geometric distance on the space of probability
distributions. This “Wasserstein” distance is thus parameterized by the
underlying ground cost between pairs of points. We refer to the monographs
[56, 49] for detailed accounts on the theory of OT, and [47] for its
computational aspects. OT distances have been used for applications as diverse
as image retrieval [48], brain imaging [22, 29], natural language processing
[34], generative models [4], and biological cell tracking [25, 50].

#### Entropic regularization.

Entropic regularization is the method of choice to scale OT to machine
learning problems. It approximates OT distance using Sinkhorn’s algorithm,
which has quadratic complexity and streams well on GPU architectures. This
line of ideas was put forward in the seminal paper by Cuturi [11], who also
emphasizes the smoothing effect, which is crucial when using Sinkhorn as a
loss function to train deep learning models. Another benefit of this
regularization is that it suffers less from the curse of dimensionality, as
proved in [20, 41]. This approach is also pivotal to scale our unsupervised
metric learning method to tackle for instance applications in genomics.

#### Metric Learning.

Metric learning is most often framed as the supervised problem of minimizing
(resp. maximizing) the distance between points in a same (resp. different)
class, and we refer to the surveys of existing approaches by [33, 5]. It is
necessary to restrict the class of distances to make the problem tractable.
The common option is arguably to consider the class of Mahalanobis distances,
which generalize the Euclidean distance and are equivalent to computing a
vectorial embedding of the data points. See for instance [60, 58, 13]. One can
apply these methods for histogram data, or use instead of Euclidean distances
more adapted discrepancies on the simplex, such as Chi-squared [44, 62] and
geodesic distances [35]. These methods however fail to capture the geometric
nature of the problem, where histograms correspond to discrete distributions
viewed as sums of localized Dirac masses.

#### OT Ground Metric Learning.

This geometry is leveraged by [12] by introducing the problem of supervised OT
ground metric learning and developing a nearest-neighbor based algorithm to
solve it. This approach is further refined in [57], which drops the triangular
inequality constraint (as we do in our approach). It is possible to restrict
the class of ground metrics, for instance using Mahalanobis [61, 31] or
geodesic distances [26] to develop more efficient learning schemes. [64]
simultaneously perform ground metric learning and matrix factorization, and
this finds applications to NLP [28]. Metric learning can also be performed
through adversarial optimization, where the metric is maximized over to
perform generative model training [19], discriminant analysis [17] and to
define robust transportation distances [45, 43]. Note that when imposing only
convex constraints, adversarial ground metric learning is a concave
maximization problem which finds applications in the modeling of crowd
congestion phenomena [6]. Another related question is the inverse problem of
estimating a ground cost from the observation of matchings or couplings [18,
53, 37, 46]. This supervised metric learning problem can be regularized using
sparsity or low-rank constraints, as explained in [14, 7].

###  1.2 Contributions

Our main contribution is the introduction in Section 3 (resp. in Section 4) of
Wasserstein eigenvectors (resp. singular vectors) as the positive eigenvectors
of a monotone homogeneous “distance map”. The associated theoretical
contributions, Theorem 1 and Proposition 2, state conditions ensuring the
existence and uniqueness of such an eigenvector.

Section 5 explains how to scale and parallelize this method by leveraging
entropic regularization through the Sinkhorn algorithm. We also show that in
the large regularization limit, our method computes metrics associated to 1-D
and 2-D embeddings along the leading principal component axes.

Numerical simulations showcase the use of these Wasserstein eigenvectors for
unsupervised metric learning on simple datasets. In Section 6 we apply our
method to a single-cell RNA sequencing dataset. We observe that the inferred
ground metric not only leads to improved clustering, but also encodes
biologically relevant information.

The Python code to reproduce the numerical results of the paper is available
online333https://github.com/gjhuizing/wasserstein-eigenvectors. The OT
computations for the results of Sections 3 and 4 (resp. 6) are computed using
the linear solver (resp. Sinkhorn’s algorithm) of the POT library [16].

##  2 Ground Costs and OT

The simplex (set of histograms) is defined as
$\Sigma_{n}\triangleq\left\\{a\in\mathbb{R}_{+}^{n}\;;\;\sum_{i}a_{i}=1\right\\}$.
The most straightforward distance between histograms is the total variation
norm $|\\!|a-b|\\!|_{1}\triangleq\sum_{i}|a_{i}-b_{i}|$. In sharp contrast,
optimal transport (Wasserstein) distances define geometrical distances by
lifting a ground cost $C$ between the bins of the histogram into a cost
between the histograms.

The ground cost $C$ should be symmetric, non-negative and with zero diagonal

$\mathcal{D}_{n}\triangleq\left\\{C\in\mathbb{R}_{+}^{n\times
n}\;;\;C^{\top}=C\quad\text{and}\quad\operatorname{diag}(C)=0_{n}\right\\}.$

The relative interior of $\mathcal{D}_{n}$ is composed of positive costs

$\mathcal{D}_{n}^{\circ}\triangleq\left\\{C\in\mathcal{D}_{n}\;;\;\forall\,i\neq
j,C_{i,j}>0\right\\}.$

The Wasserstein cost between $(a,b)\in(\Sigma_{n})^{2}$ is then obtained by
solving the linear program

$W_{C}(a,b)\triangleq\min_{P\in\Pi_{a,b}}\langle
C,\,P\rangle=\sum_{i,j}C_{i,j}P_{i,j},$ (1)

where the set of couplings between two histograms is

$\Pi_{a,b}\triangleq\left\\{P\in\mathbb{R}_{+}^{n\times
n}\;;\;P\mathds{1}_{n}=a,P^{\top}\mathds{1}_{n}=b\right\\}$

where $\mathds{1}_{n}\triangleq(1,\ldots,1)^{\top}\in\mathbb{R}^{n}$.

On $\mathcal{D}_{n}$, we compare two costs using the $\ell^{\infty}$ norm
$|\\!|C-C^{\prime}|\\!|_{\infty}$. To compare positive eigenvectors (which are
defined up to a constant) on $\mathcal{D}_{n}^{\circ}$, we use the Hilbert
projective metric

$d_{\mathcal{H}}(C,C^{\prime})\triangleq\max_{i\neq
j}(\log(C_{i,j}/C^{\prime}_{i,j}))-\min_{i\neq
j}(\log(C_{i,j}/C^{\prime}_{i,j})).$

##  3 Wasserstein Eigenvectors

We consider the problem of inferring a ground metric in an unsupervised way
from a collection of probability distributions. The rationale is that the
ground metric one uses to compare these histograms should be equal (up to an
eigenvalue constant) to the distance between them.

###  3.1 Wasserstein Distance Map

We consider a set of $m$ histograms of size $n$,
$A\triangleq(a_{j}\in\Sigma_{n})_{j=1}^{m}$. In the following, we always
assume $a_{i}\neq a_{j}$ for $i\neq j$. These histograms induce a distance map
$\Phi_{A}:\mathcal{D}_{n}\rightarrow\mathcal{D}_{m}$ defined, for
$i,j=1,\ldots,m$

$[\Phi_{A}(C)]_{i,j}\triangleq
W_{C}(a_{i},a_{j})+\tau|\\!|C|\\!|_{\infty}|\\!|a_{i}-a_{j}|\\!|_{1}.$

The parameter $\tau\geqslant 0$ is intended to play the role of a
regularization which is important for some applications.

The following proposition, whose proof can be found in Appendix A, states the
main properties of this map.

######  Proposition 1.

The map $\Phi_{A}$ is continuous from $\mathcal{D}_{n}$ to $\mathcal{D}_{m}$
and from $\mathcal{D}_{n}^{\circ}$ to $\mathcal{D}_{m}^{\circ}$. Furthermore,
it is is positively 1-homogeneous and monotone:

$\Phi_{A}(\gamma C)=\gamma\Phi_{A}(C)\quad\text{for}\quad\gamma\geqslant 0$
$C\leqslant C^{\prime}\Rightarrow\Phi_{A}(C)\leqslant\Phi_{A}(C^{\prime}).$

Lastly, it satisfies $\forall\,(C,C^{\prime})\in(\mathcal{D}_{n})^{2}$

$|\\!|\Phi_{A}(C)-\Phi_{A}(C^{\prime})|\\!|_{\infty}\leqslant(1+2\tau)|\\!|C-C^{\prime}|\\!|_{\infty}\quad\text{and}$
$\forall\,(C,C^{\prime})\in(\mathcal{D}_{n}^{\circ})^{2},\>d_{\mathcal{H}}(\Phi_{A}(C),\Phi_{A}(C^{\prime}))\leqslant
d_{\mathcal{H}}(C,C^{\prime}).$

###  3.2 Eigenvectors

The set of homogeneous monotone functions is a class of mappings which plays
an important role in topics ranging from biological modeling to stochastic
game theory and dynamic programming (Bellman and Shapley operators). We refer
to [36] for an overview.

In the case $n=m$ which we now consider, the fundamental question is to
understand the behavior of iterated compositions of such a map, which
corresponds to the power iterations method detailed in Section 3.3. This
raises the question of existence and uniqueness of eigenvectors which are
$C\in\mathcal{D}_{n}$ such that $\Phi_{A}(C)=\lambda C$ for $\lambda\geqslant
0$. Note that the $\ell^{\infty}$ bound of Proposition 1 shows that all
eigenvalues of $\Phi_{A}$ are smaller than $1+2\tau$.

The study of eigenvectors of $\Phi_{A}$ should be understood as a
generalization of the study of the eigenvectors of positive linear maps (such
as $A$) which is solved via the classical Perron-Frobenius theory [36].

#### Existence of eigenvectors.

Since $\Phi_{A}$ maps $\mathcal{D}_{n}^{\circ}$ to itself, a question of both
theoretical and practical interest is to ensure the existence of positive
(i.e. non-vanishing) eigenvectors in $\mathcal{D}_{n}^{\circ}$. The following
proposition ensures that this is the case when $\tau>0$.

######  Theorem 1.

If $\tau>0$, there exists a positive eigenvector of $\Phi_{A}$. If
$\tau\geqslant 0$, the eigenvalues associated to positive eigenvectors are all
equal.

Its proof can be found in Appendix B. Extending the existence part of this
theorem to the case $\tau=0$ is an open problem, and is out-of-reach using
classical non-linear Perron-Frobenius theorems such as [2, 1], which do not
apply. The main obstacle is that for $\tau=0$, $\Phi_{A}$ does not map
$\mathcal{D}_{n}$ to $\mathcal{D}_{n}^{\circ}$, and $0$ is actually an
eigenvalue since sparse costs are mapped to $0$ by $\Phi_{A}$.

#### Uniqueness of the eigenvector.

The second question is the uniqueness of non-negative eigenvectors. If $A$ is
too sparse and $\tau=0$, one cannot hope to have uniqueness of the leading
eigenvector. For instance, if $A$ is a permutation matrix like
$A=\mathrm{Id}_{n}$ (i.e. the histograms are Dirac masses), then all costs are
fixed points ($\Phi_{A}(C)=C$ for any $C$) so they are all leading
eigenvectors.

It does not seem obvious to provide an a priori condition (depending only on
$A$) ensuring uniqueness, but the following proposition (proved in Appendix C)
gives an a posteriori way to check the uniqueness of an eigenvector inside
$\mathcal{D}_{n}$.

######  Proposition 2.

Let $C\in\mathcal{D}_{n}$ be an eigenvector of $\Phi_{A}$. We denote
$P(a,a^{\prime})$ an optimal coupling solution of (1) for the cost $C$. If
there is a choice of optimal coupling such that the graph
$(i,j)\leftrightarrow(k,\ell)$ when $P(a_{i},a_{j})_{k,\ell}>0$ is strongly
connected, then $C$ is the unique eigenvector in $\mathcal{D}_{n}$.

In the numerical simulations showed in Section 3.5, 4.4 and 6, we checked a
posteriori that the computed eigenvectors are indeed unique.

###  3.3 Power Iterations Algorithm

The standard method to extract eigenvectors are “power iterations”
$C_{k+1}=\Phi_{A}(C_{k})/|\\!|C_{k}|\\!|_{\infty}$. In the case of linear
positive maps, Perron-Frobenius theory ensures the convergence toward the
unique positive eigenvector at a linear rate for the Hilbert metric.
Unfortunately, this result does not hold in general for the case of non-linear
maps, and $\Phi_{A}$ is only non-expansive (and not necessarily
contractacting). In practice, numerical simulations (see sections 3.5 and 4.4)
suggest that linear rates still hold for $\Phi_{A}$ when a positive singular
vector exists and is unique.

###  3.4 Examples

#### Simple cases.

While in general computing in closed form eigenvector is out of reach, it is
possible in some simple cases:

  * –

Large regularization. As $\tau\rightarrow+\infty$, then $\Phi_{A}(C)/\tau$
tends to $|\\!|C|\\!|_{\infty}T$ where $T_{i,j}=|\\!|a_{i}-a_{j}|\\!|_{1}$,
whose only eigenvector is $T$ with associated eigenvalue
$|\\!|T|\\!|_{\infty}$.

  * –

Identical histograms. If the $a_{i}$ are all equal, then $\Phi_{A}(C)=0$ so
that $0$ is the only eigenvalue.

  * –

Dirac histograms. If the $a_{i}$ are distinct Dirac masses, then $A$ is a
permutation matrix, so that $\Phi_{A}(C)=A^{\top}CA$ is a non-negative linear
map. One can apply Perron-Frobenius theory to this map, which shows that in
general it does not have a positive eigenvector. If $A=\mathrm{Id}_{n}$, then
$\Phi_{A}(C)=C$ and any $C\in\mathcal{D}_{n}$ is a leading eigenvector with
eigenvalue 1.

  * –

Pairs of histograms. When $A\in\mathbb{R}^{2\times 2}$, there is a single
eigenvector $C=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}$, and the associated
(thus dominant) eigenvalue is $|\\!|a_{1}-a_{2}|\\!|_{1}(1/2+\tau)$.

#### Block-wise interactions.

The following proposition (proved in Appendix D) shows that for a block-
diagonal histogram matrix, the cost is itself block diagonal with constant
blocks and zero on the diagonal.

######  Proposition 3.

We assume that $A=\operatorname{diag}(A_{i})$ is a block-diagonal matrix,
where $A_{i}\in\mathbb{R}_{+}^{n_{i}\times m_{i}}$. Then costs of the form
$C=(c_{i,j}\mathds{1}_{n_{i}\times m_{j}})_{i,j}$ for
$c_{i,j}\in\mathbb{R}_{+}^{*}$ and $c_{i,i}=0$ are dominant eigenvectors with
associated eigenvalue $1$.

Note that in the case of only two clusters the dominant eigenvector is
independent of the input $A$ since by symmetry, it is simply
$C=\begin{pmatrix}0_{n_{1}\times n_{1}}&1_{n_{1}\times n_{2}}\\\
1_{n_{2}\times n_{1}}&0_{n_{2}\times n_{2}}\end{pmatrix}$.

This propositon shows that when the histograms exhibit clusters by having
disjoint supports, then the eigenvector problem acts as an efficient metric
learning method, which automatically amplifies the structure of the clusters
by associating a 0 cost within each cluster, and a large cost between these
clusters. This type of behavior is important for the application to data that
contains clusters, as highlighted by the examples of MNIST handwritten digits
in section 4.4 and single-cell RNA-seq data in section 6.

#### Infinite dimensional case.

While this is not the focus of this paper, one can extend Wasserstein
eigenvectors to a possibly infinite collection of probability measures
$A\triangleq(a_{x})_{x\in\mathcal{X}}$ defined on some space $\mathcal{X}$.
For instance, for $\tau=0$, in 1-D, any cost of the form $C(x,y)=c(x-y)$ with
$c$ convex (for instance $c(s)=|s|^{p}$ for any $p\geqslant 1$) is an
eigenvector of $\Phi_{A}$ with associated eigenvalue $1$. Indeed, in such a
case, the optimal transport between $a_{x}$ and $a_{y}$ is the translation by
$x-y$ and its cost is $c(x-y)$ (see [49]).

###  3.5 Numerical Illustration: Translated Histograms

For both 1-D and 2-D settings, we generate three synthetic datasets by
translating three different templates. The datasets
$A_{1},A_{2},A_{3}\in\mathbb{R}^{n\times n}$ are thus defined by
$A_{k}\triangleq(h_{k}(x_{i}-x_{j}))_{i,j}$, where $(h_{1},h_{2},h_{3})$ are
three different (periodic) functions on the 1-D and 2-D torus (we use periodic
boundary conditions), so that the grid points $x_{i}$ are either 1-D samples
($x_{i}=i/n$) or 2-D pixel locations. We use $n=100$ points in 1-D and
$n=15\times 15=225$ pixels in 2-D.

|
---|---
Figure 1: Numerical simulations on the 1-D torus: (left) histograms whose
translations on the 1-D torus form three datasets  $A_{1}$,  $A_{2}$,  $A_{3}$
(right, top) distance maps  $c_{1}$,  $c_{2}$,  $c_{3}$ from the eigenvectors
of  $\Phi_{A_{1}}$,  $\Phi_{A_{2}}$,  $\Phi_{A_{3}}$ computed by power
iterations (right, bottom) convergence rate of the power iterations, according
to the Hilbert metric $d_{\mathcal{H}}$.

Figures 1 (for 1-D) and 2 (for 2-D) display the three templates as well as the
eigenvector metric $(C_{1},C_{2},C_{3})$, which, by translation invariance of
the problem, are of the form $(C_{k})_{i,j}=(c_{k})_{j-i}$ where
$c_{k}=(C_{k})_{0,\cdot}$ are periodic 1-D (displayed as a 1-D graph) or 2-D
(displayed as an images) functions. These results demonstrate that the learned
metrics integrate geometrical properties (symmetries, multi-modalities, etc.)
of the input datasets. For unimodal Gaussian-like distributions, the
eigenvector metric is close to $|\sin(x_{i}-x_{j})|$, but it exhibits non-
monotonic behaviors for multi-modal distributions.

Figure 2: Numerical simulations on the 2-D torus: (left, top) reshaped
histograms whose translations on the 2-D torus form three datasets  $A_{1}$,
$A_{2}$,  $A_{3}$ (left, bottom) distance maps  $c_{1}$,  $c_{2}$,  $c_{3}$
from the eigenvectors of  $\Phi_{A_{1}}$,  $\Phi_{A_{2}}$,  $\Phi_{A_{3}}$
computed by power iterations (right) convergence rate of the power iterations,
according to the Hilbert metric $d_{\mathcal{H}}$.

These two figures also report in logarithmic scales the convergence rate of
power iterations according to the Hilbert metric. This speed is always linear,
suggesting that the maps $\Phi_{A_{k}}$ are contracting and that the
eigenvector is unique (which is confirmed by running several initializations
in $\mathcal{D}_{n}^{\circ}$). The contractance rates (which is the slope of
the error curves) is highly dependent on the geometry of the templates
$h_{k}$.

We use the same datasets $A_{1}$, $A_{2}$ and $A_{3}$ (in the 1-D case) to
showcase the regularizing effect of the parameter $\tau$ (figure 3).

Figure 3: Distance maps  $c_{1}$,  $c_{2}$,  $c_{3}$ associated to the
eigenvectors of  $\Phi_{A_{1}}$,  $\Phi_{A_{2}}$,  $\Phi_{A_{3}}$ in the 1-D
case, for varying values of $\tau$.

##  4 Wasserstein Singular Vectors for Metric Learning

For many applications in learning, it makes little sense to consider the
setting $n=m$. In the general case, one rather aims for “Wasserstein singular
vectors” associated to a pair of datasets

$A=(a_{j}\in\Sigma_{n})_{j=1}^{m}\quad\text{and}\quad
B=(b_{i}\in\Sigma_{m})_{i=1}^{n}.$

The most common situation (which is also the setup associated to
dimensionality reduction methods such as PCA, ICA or NMF) is when $A$ and $B$
are the rows and the columns of the same data matrix
$U\in\mathbb{R}_{+}^{n\times m}$. This is the setup considered in the
application to genomics in Section 6. The simplest choice then is to column
and row-normalize this data matrix and define

$A=U\operatorname{diag}\left(\frac{1}{U^{\top}\mathds{1}_{n}}\right)\text{ ;
}B=U^{\top}\operatorname{diag}\left(\frac{1}{U\mathds{1}_{m}}\right).$ (2)

One can also apply Sinkhorn’s matrix scaling (see Section 5.1) to define a
symmetric normalization to a bistochastic matrix.

###  4.1 Wasserstein Singular Vectors

These two sets thus define maps $\Phi_{A}$ and $\Phi_{B}$ between
$\mathcal{D}_{n}$ and $\mathcal{D}_{m}$, so that one can apply the
construction of Wasserstein eigenvectors to the maps
$\Phi_{A}\circ\Phi_{B}:\mathcal{D}_{m}\mapsto\mathcal{D}_{m}$ and
$\Phi_{B}\circ\Phi_{A}:\mathcal{D}_{n}\mapsto\mathcal{D}_{n}$, to define a
pair $(C,D)\in\mathcal{D}_{n}\times\mathcal{D}_{m}$ of Wasserstein singular
vectors, with associated singular values $(\lambda,\mu)$ such that

$D=\lambda\Phi_{A}(C)\quad\text{and}\quad C=\mu\Phi_{B}(D).$

Thanks to Proposition 1, one has $\lambda,\mu\leqslant 1+2\tau$. Note that in
the linear case, one has $\lambda=\mu$, but this is not necessarily true in
our non-linear setting. The proof of Theorem 1 carries over to singular
vectors, which ensures existence of positive singular vectors whenever
$\tau>0$.

This construction thus operates as a joint metric learning method. For
instance, in the case where $A$ and $B$ are computed from the rows and columns
of a common data matrix, singular vectors $(C,D)$ define two intrinsic metrics
between rows and columns respectively.

Proposition 3 extends to this case, when $A$ and $B$ are diagonal block
matrices (possibly with non square blocks) with the same block structures. On
these simple models, this reflects the fact that these intrinsic metric
capture the block-wise clustering structure of the data. This explains in part
the efficiency of this approach for the application in genomics showcased in
Section 6.

###  4.2 Power Iterations Algorithm

Starting from $C_{0}$, the iterations

$D_{k+1}\triangleq\frac{\Phi_{A}(C_{k})}{|\\!|\Phi_{A}(C_{k})|\\!|_{\infty}}\;\text{and}\;C_{k+1}\triangleq\frac{\Phi_{B}(D_{k+1})}{|\\!|\Phi_{B}(D_{k+1})|\\!|_{\infty}}$
(3)

are expected to converge to a pair of leading singular vectors.

###  4.3 Numerical Illustration: Mean/Scale Families

We generalize the setting of Section 3.5 by considering both mean and scale
parameters $(\mu,\sigma)$. The histogram matrices $A$ and $B$ are obtained
using the normalization (2) of the data matrix $U$ defined as
$U_{i,j}\triangleq h\left(\frac{x_{i}-\mu_{j}}{\sigma_{j}}\right)$ where
$x_{i}=i/n$ is a 1-D grid while
$(\mu_{j},\sigma_{i})\in[0,1]\times\mathbb{R}_{+}^{*}$ is a 2-D regular grid
of $m$ points over the mean/scale plane. Here $h$ is a 1-D Gaussian-like
profile on the torus (we use periodic boundary conditions). The histograms
$a_{i}$ thus form a 2-D parametric family of histograms on a 1-D domain, which
highlights the need to move from eigenvectors to singular vectors to cope with
this mismatch between the two spaces.

Figure 4 displays the singular vectors computed numerically using power
iterations (3). The first singular vector $C\in\mathcal{D}_{n}$ defines a 1-D
metric, which is by translation invariance of the problem of the form
$C_{i,i^{\prime}}=c_{i^{\prime}-i}$. The left panel of the figure thus shows
the 1-D function $c=C_{0,\cdot}$ (which we found numerically to be
approximately equal to $c_{i}=|\sin(x_{i})|$). The second singular vector
$D\in\mathcal{D}_{m}$ defines a metric on the 2-D mean/scale plane, and we
display the associated distance maps $(D_{j_{0},\cdot},D_{j_{1},\cdot})$ for
two different indexes (corresponding to zero mean $\mu_{j_{1}}=\mu_{j_{2}=0}$
bumps with respectively a small width $\sigma_{j_{1}}$ and a large width
$\sigma_{j_{2}}$). These distance maps are vizualized as 2-D images over this
plane, which highlights how the fixed point iteration is able to integrate the
information between the whole family to define a meaningful notion of cost
between Gaussian-like functions. Note that on $\mathbb{R}$, the Wasserstein
distance for the cost $|x-x^{\prime}|^{2}$ between two Gaussians of parameters
$(\mu,\sigma)$ and $(\mu^{\prime},\sigma^{\prime})$ is known to be the
Euclidean distance on the mean/scale plane
$|\\!|(\mu,\sigma)-(\mu^{\prime},\sigma^{\prime})|\\!|^{2}$ [54], but this
numerical simulation indicates that the fixed point metric behaves quite
differently.

|  |
---|---|---
$C_{0,\cdot}$ | $D_{j_{0},\cdot}$ | $D_{j_{1},\cdot}$
Figure 4: Representation of the singular vectors $(C,D)$ associated to
translated and scaled bump functions.

###  4.4 Numerical Illustration: MNIST

We consider a small subset of $m=55$ samples of the MNIST handwritten digits
dataset, which are gray-scale images images sampled at $n=28^{2}=784$ pixels
locations $(x_{i})_{i=1}^{n}$. We select only two classes (zeros and ones) in
order to interpret the singular vectors more easily. We can thus represent the
dataset as a matrix $U\in\mathbb{R}_{+}^{n\times m}$. The histogram matrices
$A$ and $B$ are obtained using the normalization (2).

The operators $\Phi_{A}$ and $\Phi_{B}$ are computed for $\tau=0$. The
singular vectors $(C\in\mathcal{D}_{n},D\in\mathcal{D}_{m})$ are computed
using the fixed point iterations power iterations (3), initialized with a
random cost matrix in $\mathcal{D}_{n}^{\circ}$.

Figure 5: Top: Euclidean distance matrix between digits, and the associated
MDS projection ; Bottom: singular vector $D$ learned by power iterations, and
the associated MDS projection.

Figure 5, bottom row, shows the pairwise distances $D_{j,j^{\prime}}$ between
all pairs of MNIST images. By the singular vector property, these distances
are Wasserstein distances for the ground cost $C_{i,i^{\prime}}$ between pairs
of pixels $x_{i}$ and $x_{i^{\prime}}$. This figure highlights how this
learned ground distance matrix improves over the usual squared Euclidean
distance $|\\!|x_{i}-x_{i^{\prime}}|\\!|^{2}$ between the pixels (shown on top
row). The classes are better separated, which is also apparent on the 2-D
embedding computed using classical multi-dimensional scaling (MDS).

Figure 6 shows ground metric distance maps $D_{i_{k},\cdot}$ for three
different reference pixels $(i_{1},i_{2},i_{3})$. These maps are visualized as
2-D images, so that $D_{i,i^{\prime}}$ is the distance between the reference
pixel $x_{i}$ and another pixel $x_{i^{\prime}}$. Note that these distance
maps are not defined (set to 0) outside the support of the digits. From left
to right, the dark pixels are: for $C_{i_{1},\cdot}$ pixels strongly
associated to digit zero ; for $C_{i_{2},\cdot}$ pixels strongly associated to
digit one ; for $C_{i_{3},\cdot}$ pixels associated equally to both classes.

Figure 6: Distance map $(C_{i_{1},\cdot},C_{i_{2},\cdot},C_{i_{3},\cdot})$ of
the ground cost $C$, for different pixels $(i_{1},i_{2},i_{3})$ (represented
using a red star).

##  5 Entropic Regularization

In order to be able to efficiently compute the $\Phi_{A}$ mapping, one can use
entropic regularization to define a meaningful smooth approximation of
$\Phi_{A}^{\varepsilon}$. All the entries of $\Phi_{A}^{\varepsilon}$ can be
computed in parallel using the Sinkhorn algorithm, which is crucial to scale
the method to the application in genomics considered in Section 6. We also
show in Section 5.2 that the limit of large entropic regularization operates
metric learning using principal component analysis (PCA) dimensionality
reduction.

###  5.1 Sinkhorn Divergences

#### Entropic approximation.

The Sinkhorn cost is defined as

$W_{C}^{\varepsilon}(a,b)\triangleq\min_{P\in\Pi_{a,b}}\langle
C,\,P\rangle+\varepsilon|\\!|C|\\!|_{\infty}H(P),$

where the negative entropy is $H(P)\triangleq\sum_{i,j}P_{i,j}\log(P_{i,j})$.
One can show that
$|W_{C}^{\varepsilon}-W_{C}|\sim\varepsilon\log(\varepsilon)$ [8]. This
approximation error can be further reduced to $\sim\varepsilon^{2}$ [10] by
using the debiased Sinkhorn divergence introduced in [19]

$\bar{W}_{C}^{\varepsilon}(a,b)\triangleq
W_{C}^{\varepsilon}(a,b)-\frac{1}{2}W_{C}^{\varepsilon}(a,a)-\frac{1}{2}W_{C}^{\varepsilon}(b,b).$

This debiasing is also crucial to ensure that
$\bar{W}_{C}^{\varepsilon}(a,a)=0$.

Similarly to the un-regularized case (corresponding to $\varepsilon=0$), we
define the distance map

$\Phi_{A}^{\varepsilon}(C)\triangleq(\bar{W}_{C}^{\varepsilon}(a_{i},a_{j})+\tau|\\!|C|\\!|_{\infty}|\\!|a_{i}-a_{j}|\\!|_{1})_{i,j}.$

Note that this map is still 1-homogeneous and it maps zero diagonal matrices
to zero diagonal matrices. It has been showed in [15] that
$C\in\mathcal{D}_{n}$ is mapped to
$\Phi_{A}^{\varepsilon}(C)\in\mathcal{D}_{m}$ under the condition that
$K\triangleq e^{-\frac{C}{|\\!|C|\\!|_{\infty}\varepsilon}}$ is a valid
positive kernel (i.e. has positive eigenvalues). While it is not clear that
such a condition is maintained during power iterations, we observed
numerically that it is still the case in practice.

#### Sinkhorn’s algorithm.

The Sinkhorn cost can be computed by the dual formula

$W_{C}^{\varepsilon}(a,b)=\varepsilon\left(\langle\log(u),\,a\rangle+\langle\log(v),\,b\rangle-\langle
Kv,\,u\rangle\right),$

where $(u,v)$ are obtained by iterating the following Sinkhorn fixed point

$u\leftarrow\frac{a}{K^{\top}v}\quad\text{and}\quad v\leftarrow\frac{b}{Ku}.$

This allows one to compute with precision $\varepsilon$ the $m^{2}$ entries of
$\Phi_{A}^{\varepsilon}(C)$ in $O((mn)^{2}/\varepsilon^{2})$ operations [3],
using a parallelizable algorithm that is well suited for GPU computations.

###  5.2 Connexion with PCA and Mahalanobis Metric Learning

For simplicity, we consider the case $\tau=0$ in this section. Using [19], one
has that the large $\varepsilon$-limit corresponds to the use of the
Mahalanobis distance associated to the kernel $-C$, which is often referred to
as Maximum Mean Discrepancy (MMD) distances [23].

######  Proposition 4.

When $\varepsilon\rightarrow+\infty$,

$\Phi_{A}^{\varepsilon}(C)_{i,j}\rightarrow\Phi_{A}^{\infty}(C)_{i,j}\triangleq\frac{1}{2}\langle-C(a_{i}-a_{j}),\,a_{i}-a_{j}\rangle.$

The map $\Phi_{A}^{\infty}$ is linear, but in sharp contrast to $\Phi_{A}$ it
is not positive in the classical sense, i.e. it is not mapping to
$\mathcal{D}_{n}$. The correct notion of positivity (or rather negativity) to
be used here are conditionally negative kernels [27], which is defined by

$C\preceq^{\star}0\quad\Leftrightarrow\quad\forall z\in U_{n},\langle
Cz,\,z\rangle\leqslant 0$

where $U_{n}\triangleq\left\\{z\in\mathbb{R}^{n}\;;\;\langle
z,\,\mathds{1}_{n}\rangle\right\\}=0$. The following proposition constructs a
sub-cone of $\mathcal{D}_{n}$ suited for the study of $\Phi_{A}^{\infty}$.

######  Proposition 5.

We define the convex cone

$\mathcal{K}_{n}\triangleq\left\\{C\in\mathbb{R}_{+}^{n\times
n}\;;\;C\preceq^{\star}0,\operatorname{diag}(C)=0_{n}\right\\}\subset\mathcal{D}_{n}^{\circ}.$

One has (see [27, Prop.2.24]) that

$C\in\mathcal{K}_{n}\quad\Leftrightarrow\quad\exists(u_{i}\in\mathbb{R}^{d})_{i=1}^{n},\>C_{i,j}=|\\!|u_{i}-u_{j}|\\!|^{2}$

where $d=d(C)\in\\{1,\ldots,n\\}$ (it depends on the matrix) is the
dimensionality the embedding. One has that $\Phi_{A}^{\varepsilon}$ maps
$\mathcal{K}_{n}$ to itself.

This cone is thus precisely the cone of Euclidean distance matrices. Since for
$C\neq 0$, $\Phi_{A}(C)\neq 0$ (0 is not an eigenvalue), this ensure that
there exists eigenvectors of $\Phi_{A}^{\infty}$ in $\mathcal{K}_{n}$ (which
are thus positive in the classical sense). Unfortunately, the linear map
$\Phi_{A}^{\varepsilon}$ does not map $\mathcal{K}_{n}$ to its interior, so
that one cannot apply Perron-Frobenius theory to ensure uniqueness of
eigenvectors in $\mathcal{K}_{n}$. The following proposition, proved in
Appendix E, indeed shows that one can construct several distinct eigenvectors.
They correspond to the computation of squared Euclidean costs $C$ for “PCA”
projections of the data in dimension $d(C)=1$ or $d(C)=2$.

######  Proposition 6.

We denote $\tilde{a}_{i}\triangleq a_{i}-\frac{1}{n}\sum_{j}a_{j}$ and
$\tilde{A}$ the associated centered dataset. For every eigenvector
$u\in\mathbb{C}^{n}$ of $\tilde{A}^{\top}$ with eigenvalue
$\lambda\in\mathbb{C}$, the matrix
$C=(|u_{i}-u_{j}|^{2})_{i,j}\in\mathcal{K}_{n}$ ($|\cdot|$ being the modulus
of complex numbers) is an eigenvector of $\Phi_{A}^{\infty}$ with eigenvalue
$2|\lambda(u)|^{2}$.

This proposition thus constructs non-negative eigenvectors $C$ of
$\Phi_{A}^{\infty}$ which correspond to squared Euclidean distance matrices in
a 1-D (resp. 2-D) embedding space when $u$ is real (resp. complex). This
analysis extends to the singular vector setting exposed in Section 4, and in
this case these embeddings are exactly PCA coordinates. Note however that this
analysis does not rule out the existence of other positive eigenvectors
associated to higher dimensional embeddings.

##  6 Metric Learning for Single Cell Genomics

#### Context.

RNA sequencing (RNA-seq) is a high-throughput sequencing technology enabling
the measurement of gene expression levels in hundreds of biological samples
[51]. The analysis of RNA-seq data has improved our understanding of cancer
heterogeneity and has had a big clinical impact in the context of precision
medicine [59]. RNA-seq data corresponds to the average gene expression of all
the cells present in a biological sample, but recent technological advances
have enabled the sequencing of RNA on individual cells [52]. The analysis of
single-cell RNA-seq (scRNA-seq) data has offered unprecedented insights in
cellular heterogeneity and disease mechanisms [55, 63].

A scRNA-seq dataset can be represented as a sparse matrix of expression levels
with genes on rows and cells on columns. Classical “bin-to-bin” distances
between cells (e.g. the Euclidean loss) do not carry any biological
information beyond differences in gene expression values across cells. This
makes these distances susceptible to noise and highly dependant on data
preprocessing [32]. In addition, the phenotype of a cell is determined by its
underlying gene-gene relationships (regulatory mechanisms) which are ignored
by bin-to-bin distances. In contrast, our unsupervised ground metric learning
method offers a way to take these relationships into account.

Figure 7: Top: Euclidean distance matrix between cells and the associated MDS
projection ; Bottom: singular vector distance matrix $D$ between cells, and
the associated MDS projection.

#### Dataset.

We use scRNA-seq data from [39], which contains 206 cells from three different
cell lines, i.e. three different populations of cells. We log-transform the
integer counts ($x\to\log(1+x)$) and select the 1000 most varying genes. This
preprocessing reflects common practices outlined in [40]. The preprocessed
dataset can be represented as a matrix $U\in\mathbb{R}_{+}^{n\times m}$ where
$n=206$ is the number of cells and $m=1000$ is the number of genes. The
histogram matrices $A$ and $B$ are obtained through the normalization (2).

#### Numerical experiments.

We compute the singular vectors $C$ and $D$ using the power iterations
described earlier. The pairwise Sinkhorn divergences are computed using POT
[16], with regularization parameters $\varepsilon=10^{-2}$ and $\tau=0$. We
compare the resulting distance matrix between cells with the Euclidean
distance matrix and analyse the gene-to-gene ground metric. The idea here is
to evaluate the Euclidean and OT distance matrices based on their capability
to capture the similarity between cells extracted from the same cell line and
the separation across cells extracted from different cell lines.

#### Results.

Figure 7 shows the Euclidean distance matrix between cells, the singular
vector $D$, and the corresponding MDS projections. We observe that the OT
distance matrix with the learned ground metric better separates the different
cell lines. For visualization purposes, we apply agglomerative hierarchical
clustering on the learned cost matrix $C$ between genes. Sorting the genes by
cluster allows us to observe clear similarities in the cost matrix (figure 8).
These similarities correspond to the similar patterns of expression across the
various cell lines. Figure 9 provides a visualization of the dataset where
these patterns are clearly visible.

Intuitively, genes with similar expression patterns are treated by OT
distances as interchangeable. It follows that to be clustered together, two
cells do not need to express the exact same genes, but rather genes with
similar behaviour. This motivates the use of optimal transport for
unsupervised clustering tasks on noisy data like scRNA-seq data.

Figure 8: Singular vector cost matrix $C$ between genes. For visualization
purposes, genes are sorted by hierarchical clustering. Figure 9: Dataset, with
genes arranged according to clustering of singular vector $C$

##  7 Conclusion and Perspectives

In this paper we defined Wasserstein eigenvectors and provided results on
their existence and uniqueness. We then extended these results to singular
vectors, which are well suited to machine learning problems. We showed that
unsupervised ground metric learning could be achieved by computing such
singular vectors. This method is able to scale to large datasets using the
highly parallelizable Sinkhorn algorithm. In the large regularization limit,
our method is linked to principal component analysis. We applied our method to
synthetic data, to the MNIST handwritten digits dataset, and to a single-cell
RNA sequencing dataset. In all cases, the ground metric learned iteratively is
intuitively interpretable. In particular, the ground metric learned on
biological data not only leads to improved clustering, but also encodes
biologically relevant information.

Theoretical perspectives include further results on the existence of positive
eigenvectors, in particular for $\tau=0$ and for $\varepsilon>0$. In addition,
integrating unbalanced optimal transport [38, 9] into the method could avoid
the need for the step of normalization to histograms. Applying our method to
large single cell datasets is also a promising avenue to extend the
applicability of OT to new classes of problems in genomics.

## Acknowledgements

We thank Stéphane Gaubert for very useful advises on non-linear Perron-
Frobenius theory. The project leading to this publication has received funding
from the Agence Nationale de la Recherche (ANR) project scMOmix and Sanofi
iTech Awards. The work of G. Peyré is supported by the European Research
Council (ERC project NORIA) and by the French government under management of
Agence Nationale de la Recherche as part of the “Investissements d’avenir”
program, reference ANR19-P3IA-0001 (PRAIRIE 3IA Institute).

## References

  * [1] Marianne Akian, Stéphane Gaubert, and Antoine Hochart.  A game theory approach to the existence and uniqueness of nonlinear perron-frobenius eigenvectors.  arXiv preprint arXiv:1812.09871, 2018.
  * [2] Marianne Akian, Stéphane Gaubert, and Roger Nussbaum.  Uniqueness of the fixed point of nonexpansive semidifferentiable maps.  Transactions of the American Mathematical Society, 368(2):1271–1320, 2016.
  * [3] Jason Altschuler, Jonathan Weed, and Philippe Rigollet.  Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration.  In Advances in Neural Information Processing Systems, pages 1961–1971, 2017.
  * [4] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In International conference on machine learning, pages 214–223. PMLR, 2017.
  * [5] Aurélien Bellet, Amaury Habrard, and Marc Sebban.  A survey on metric learning for feature vectors and structured data.  arXiv preprint arXiv:1306.6709, 2013.
  * [6] Fethallah Benmansour, Guillaume Carlier, Gabriel Peyré, and Filippo Santambrogio.  Derivatives with respect to metrics and applications: subgradient marching algorithm.  Numerische Mathematik, 116(3):357–381, 2010.
  * [7] Guillaume Carlier, Arnaud Dupuy, Alfred Galichon, and Yifei Sun.  Sista: learning optimal transport costs under sparsity constraints.  arXiv preprint arXiv:2009.08564, 2020.
  * [8] Guillaume Carlier, Vincent Duval, Gabriel Peyré, and Bernhard Schmitzer.  Convergence of entropic schemes for optimal transport and gradient flows.  SIAM Journal on Mathematical Analysis, 49(2):1385–1418, 2017.
  * [9] Lénaïc Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  Unbalanced optimal transport: Dynamic and kantorovich formulations.  Journal of Functional Analysis, 274(11):3090–3123, 2018.
  * [10] Lenaic Chizat, Pierre Roussillon, Flavien Léger, François-Xavier Vialard, and Gabriel Peyré.  Faster wasserstein distance estimation with the sinkhorn divergence.  In Proc. NeurIPS’20, 2020.
  * [11] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In Adv. in Neural Information Processing Systems, pages 2292–2300, 2013.
  * [12] Marco Cuturi and David Avis.  Ground metric learning.  The Journal of Machine Learning Research, 15(1):533–564, 2014.
  * [13] Jason V Davis and Inderjit S Dhillon.  Structured metric learning for high dimensional problems.  In Proceedings of the 14th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 195–203, 2008.
  * [14] Arnaud Dupuy, Alfred Galichon, and Yifei Sun.  Estimating matching affinity matrices under low-rank constraints.  Information and Inference: A Journal of the IMA, 8(4):677–689, 2019\.
  * [15] Jean Feydy, Thibault Séjourné, François-Xavier Vialard, Shun-ichi Amari, Alain Trouvé, and Gabriel Peyré.  Interpolating between optimal transport and mmd using sinkhorn divergences.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2681–2690, 2019.
  * [16] Rémi Flamary and Nicolas Courty.  Pot python optimal transport library, 2017.
  * [17] Rémi Flamary, Marco Cuturi, Nicolas Courty, and Alain Rakotomamonjy.  Wasserstein discriminant analysis.  Machine Learning, 107(12):1923–1945, 2018.
  * [18] Alfred Galichon and Bernard Salanié.  Cupid’s invisible hand: Social surplus and identification in matching models.  Available at SSRN 1804623, 2020.
  * [19] A. Genevay, G. Peyré, and M. Cuturi.  Learning generative models with sinkhorn divergences.  In Proc. AISTATS’18, pages 1608–1617, 2018.
  * [20] Aude Genevay, Lénaic Chizat, Francis Bach, Marco Cuturi, and Gabriel Peyré.  Sample complexity of sinkhorn divergences.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1574–1583. PMLR, 2019.
  * [21] Alison L Gibbs and Francis Edward Su.  On choosing and bounding probability metrics.  International statistical review, 70(3):419–435, 2002.
  * [22] Alexandre Gramfort, Gabriel Peyré, and Marco Cuturi.  Fast optimal transport averaging of neuroimaging data.  In International Conference on Information Processing in Medical Imaging, pages 261–272. Springer, 2015.
  * [23] Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  Journal of Machine Learning Research, 13(Mar):723–773, 2012.
  * [24] Jeremy Gunawardena and Mike Keane.  On the existence of cycle times for some nonexpansive maps.  Technical report, Citeseer, 1995.
  * [25] Tatsunori Hashimoto, David Gifford, and Tommi Jaakkola.  Learning population-level diffusions with generative rnns.  In International Conference on Machine Learning, pages 2417–2426. PMLR, 2016.
  * [26] Matthieu Heitz, Nicolas Bonneel, David Coeurjolly, Marco Cuturi, and Gabriel Peyré.  Ground metric learning on graphs.  Journal of Mathematical Imaging and Vision, pages 1–19, 2020.
  * [27] Thomas Hofmann, Bernhard Schölkopf, and Alexander J Smola.  Kernel methods in machine learning.  The annals of statistics, pages 1171–1220, 2008.
  * [28] Gao Huang, Chuan Quo, Matt J Kusner, Yu Sun, Kilian Q Weinberger, and Fei Sha.  Supervised word mover’s distance.  In Proceedings of the 30th International Conference on Neural Information Processing Systems, pages 4869–4877, 2016.
  * [29] Hicham Janati, Marco Cuturi, and Alexandre Gramfort.  Spatio-temporal alignments: Optimal transport through space and time.  In International Conference on Artificial Intelligence and Statistics, pages 1695–1704. PMLR, 2020.
  * [30] Leonid Kantorovich.  On the transfer of masses (in Russian).  Doklady Akademii Nauk, 37(2):227–229, 1942.
  * [31] Tanguy Kerdoncuff, Rémi Emonet, and Marc Sebban.  Metric Learning in Optimal Transport for Domain Adaptation.  In International Joint Conference on Artificial Intelligence, Kyoto, Japan, January 2021.
  * [32] Taiyun Kim, Irene Rui Chen, Yingxin Lin, Andy Yi-Yang Wang, Jean Yee Hwa Yang, and Pengyi Yang.  Impact of similarity metrics on single-cell rna-seq data clustering.  Briefings in bioinformatics, 20(6):2316–2326, 2019.
  * [33] Brian Kulis et al.  Metric learning: A survey.  Foundations and trends in machine learning, 5(4):287–364, 2012\.
  * [34] Matt Kusner, Yu Sun, Nicholas Kolkin, and Kilian Q Weinberger.  From word embeddings to document distances.  In Proc. of the 32nd Intern. Conf. on Machine Learning, pages 957–966, 2015.
  * [35] Tam Le and Marco Cuturi.  Unsupervised riemannian metric learning for histograms using aitchison transformations.  In International Conference on Machine Learning, pages 2002–2011. PMLR, 2015.
  * [36] Bas Lemmens and Roger Nussbaum.  Nonlinear Perron-Frobenius Theory, volume 189.  Cambridge University Press, 2012.
  * [37] Ruilin Li, Xiaojing Ye, Haomin Zhou, and Hongyuan Zha.  Learning to match via inverse optimal transport.  Journal of machine learning research, 20, 2019.
  * [38] Matthias Liero, Alexander Mielke, and Giuseppe Savaré.  Optimal entropy-transport problems and a new hellinger–kantorovich distance between positive measures.  Inventiones mathematicae, pages 1–149, 2015.
  * [39] Longqi Liu, Chuanyu Liu, Andrés Quintero, Liang Wu, Yue Yuan, Mingyue Wang, Mengnan Cheng, Lizhi Leng, Liqin Xu, Guoyi Dong, et al.  Deconvolution of single-cell multi-omics layers reveals regulatory heterogeneity.  Nature communications, 10(1):1–10, 2019.
  * [40] Malte D Luecken and Fabian J Theis.  Current best practices in single-cell rna-seq analysis: a tutorial.  Molecular systems biology, 15(6):e8746, 2019.
  * [41] Gonzalo Mena and Jonathan Weed.  Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem.  arXiv preprint arXiv:1905.11882, 2019.
  * [42] G. Monge.  Mémoire sur la théorie des déblais et des remblais.  Histoire de l’Académie Royale des Sciences, pages 666–704, 1781\.
  * [43] Jonathan Niles-Weed and Philippe Rigollet.  Estimation of wasserstein distances in the spiked transport model.  arXiv preprint arXiv:1909.07513, 2019.
  * [44] Samyeul Noh.  $\chi^{2}$ metric learning for nearest neighbor classification and its analysis.  In Proceedings of the 21st International Conference on Pattern Recognition (ICPR2012), pages 991–995. IEEE, 2012.
  * [45] François-Pierre Paty and Marco Cuturi.  Subspace robust wasserstein distances.  In International Conference on Machine Learning, pages 5072–5081. PMLR, 2019.
  * [46] François-Pierre Paty and Marco Cuturi.  Regularized optimal transport is ground cost adversarial.  In International Conference on Machine Learning, pages 7532–7542. PMLR, 2020.
  * [47] Gabriel Peyré, Marco Cuturi, et al.  Computational optimal transport.  Foundations and Trends® in Machine Learning, 11(5-6):355–607, 2019.
  * [48] 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, November 2000.
  * [49] Filippo Santambrogio.  Optimal Transport for applied mathematicians, volume 87 of Progress in Nonlinear Differential Equations and their applications.  Springer, 2015.
  * [50] Geoffrey Schiebinger, Jian Shu, Marcin Tabaka, Brian Cleary, Vidya Subramanian, Aryeh Solomon, Joshua Gould, Siyan Liu, Stacie Lin, Peter Berube, et al.  Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming.  Cell, 176(4):928–943, 2019.
  * [51] Rory Stark, Marta Grzelak, and James Hadfield.  Rna sequencing: the teenage years.  Nature Reviews Genetics, 20(11):631–656, 2019.
  * [52] Oliver Stegle, Sarah A Teichmann, and John C Marioni.  Computational and analytical challenges in single-cell transcriptomics.  Nature Reviews Genetics, 16(3):133–145, 2015.
  * [53] Andrew M Stuart and Marie-Therese Wolfram.  Inverse optimal transport.  SIAM Journal on Applied Mathematics, 80(1):599–619, 2020.
  * [54] Asuka Takatsu et al.  Wasserstein geometry of gaussian measures.  Osaka Journal of Mathematics, 48(4):1005–1026, 2011.
  * [55] Amos Tanay and Aviv Regev.  Scaling single-cell genomics from phenomenology to mechanism.  Nature, 541(7637):331–338, 2017.
  * [56] Cedric Villani.  Topics in C. Transportation.  Graduate studies in Math. AMS, 2003.
  * [57] Fan Wang and Leonidas J Guibas.  Supervised earth mover’s distance learning and its computer vision applications.  In European Conference on Computer Vision, pages 442–455. Springer, 2012.
  * [58] Kilian Q Weinberger, John Blitzer, and Lawrence K Saul.  Distance metric learning for large margin nearest neighbor classification.  In Advances in neural information processing systems, pages 1473–1480, 2006.
  * [59] John N Weinstein, Eric A Collisson, Gordon B Mills, Kenna R Mills Shaw, Brad A Ozenberger, Kyle Ellrott, Ilya Shmulevich, Chris Sander, and Joshua M Stuart.  The cancer genome atlas pan-cancer analysis project.  Nature genetics, 45(10):1113–1120, 2013.
  * [60] Eric P Xing, Andrew Y Ng, Michael I Jordan, and Stuart Russell.  Distance metric learning with application to clustering with side-information.  In NIPS, volume 15, page 12, 2002.
  * [61] Jie Xu, Lei Luo, Cheng Deng, and Heng Huang.  Multi-level metric learning via smoothed wasserstein distance.  In IJCAI, pages 2919–2925, 2018.
  * [62] Wei Yang, Luhui Xu, Xiaopan Chen, Fengbin Zheng, and Yang Liu.  Chi-squared distance metric learning for histogram data.  Mathematical Problems in Engineering, 2015, 2015.
  * [63] Guo-Cheng Yuan, Long Cai, Michael Elowitz, Tariq Enver, Guoping Fan, Guoji Guo, Rafael Irizarry, Peter Kharchenko, Junhyong Kim, Stuart Orkin, et al.  Challenges and emerging directions in single-cell analysis.  Genome biology, 18(1):1–8, 2017.
  * [64] Gloria Zen, Elisa Ricci, and Nicu Sebe.  Simultaneous ground metric learning and matrix factorization with earth mover’s distance.  In 2014 22nd International Conference on Pattern Recognition, pages 3690–3695. IEEE, 2014.

##  Appendix A Proof of Proposition 1

The homogeneity and monotony of $\Phi_{A}$ follow from its definition. Note
that $\Phi_{A}$ is a vector-valued concave function (each coordinate being an
infimum of linear forms) and hence it is continuous on $\mathbb{R}^{n\times
n}$ (and actually Lipschitz for $\ell^{\infty}$ and Hilbert metrics as we now
show). Non-expansiveness for the Hilbert metric is true for any such map, see
for instance [24]. Lipschitz continuity for the $\ell^{\infty}$ norm follows
from

$\displaystyle|W_{C}(a,a^{\prime})-W_{C^{\prime}}(a,a^{\prime})|$
$\displaystyle=|\min_{\Pi_{a,a^{\prime}}}\langle
C,\,P\rangle-\min_{\Pi_{a,a^{\prime}}}\langle C^{\prime},\,P^{\prime}\rangle|$
$\displaystyle\leqslant\max_{P\in\Pi_{a,a^{\prime}}}\langle|C-C^{\prime}|,\,P\rangle$
$\displaystyle\leqslant|\\!|C-C^{\prime}|\\!|_{\infty}|\\!|P|\\!|_{1}$
$\displaystyle=|\\!|C-C^{\prime}|\\!|_{\infty},$

since $|\min(u)-\min(v)|\leqslant\max|u-v|$ and $|\\!|P|\\!|_{1}=1$.

##  Appendix B Proof of Theorem 1

#### Proof of existence of eigenvectors.

We consider

$\bar{\Phi}_{A}(C)_{i,j}\triangleq\frac{\Phi_{A}(C)_{i,j}}{|\\!|\Phi_{A}(C)|\\!|_{\infty}}$

which is well defined and continuous when $\Phi_{A}(C)\neq 0$. For $\rho>0$,
we consider the compact set

$\mathcal{K}_{\rho}\triangleq\left\\{C\in\mathcal{D}_{n}\;;\;|\\!|C|\\!|_{\infty}=1,\forall\,i\neq
j,C_{i,j}\geqslant\rho\right\\}.$

A classical result (see for instance [21]) states that for
$C\in\mathcal{K}_{\rho}$

$\frac{\rho}{2}|\\!|a-b|\\!|_{1}\leqslant
W_{C}(a,b)\leqslant\frac{1}{2}|\\!|a-b|\\!|_{1}.$

Denoting

$\gamma\triangleq\frac{\min_{i\neq
j}|\\!|a_{i}-a_{j}|\\!|_{1}}{\max_{i,j}|\\!|a_{i}-a_{j}|\\!|_{1}}\leqslant 1$

one thus has for $i\neq j$

$\bar{\Phi}_{A}(C)_{i,j}\geqslant\gamma\frac{\rho/2+\tau}{1/2+\tau}.$

This shows that for

$0<\rho\leqslant\frac{2\gamma\tau}{2\tau+1-\gamma}$

one has $\bar{\Phi}_{A}(C)\in\mathcal{K}_{\rho}$. So for such $\rho$,
$\bar{\Phi}_{A}$ is a continuous map from the locally contractible compact
$\mathcal{K}_{\rho}$ to itself, so using Brouwer theorem, it has a fixed
point, which is an eigenvector of $\Phi_{A}$.

#### Proof of equality of the eigenvalues.

We actually show a stronger statement, namely that if $\lambda>0$ is such that
there exists $C\in\mathcal{D}_{n}^{\circ}$ with $\Phi_{A}(C)=\lambda C$, then
it satisfies for any $C^{\prime}$ and $i\neq j$

$\lim_{k\rightarrow+\infty}\frac{1}{k}\log(\Phi^{k}(C^{\prime})_{i,j})=\log(\lambda),$

where $\Phi_{A}^{k}=\Phi_{A}\circ\ldots\circ\Phi_{A}$ (composed $k$-times). In
particular, such a $\lambda\in]0,1]$ is unique. Indeed, using the non-
expansivity of Proposition 1 (which is also true for the Thompson metric), one
has

$|\\!|k\log(\lambda)+\log(C)-\log(\Phi^{k}_{A}(C^{\prime}))|\\!|_{\infty}\leqslant|\\!|\log(C/C^{\prime})|\\!|_{\infty}$

and hence the result after division by $k$.

##  Appendix C Proof of Proposition 2

We make use of Theorem 7.5 of [2]. It requires that the semi-differential of
$\Phi_{A}$ at $C$ has itself a unique positive eigenvector in
$\mathcal{D}_{n}$ (and in fact this eigenvector is also $C$). From the
enveloppe theorem, upper-gradients of the concave function $C\mapsto
W_{C}(a,a^{\prime})$ are the elements $P(a,a^{\prime})$ (so if the optimal
coupling is unique, then this map is differentiable). The semi-differential of
$\Phi_{A}$ at $C$ is thus composed by the of positive linear operators of the
form

$\Psi_{C}:\Delta\in\mathcal{D}_{n}\mapsto(\langle\Delta,\,P(a_{i},a_{j})\rangle)_{i,j}\in\mathcal{D}_{n}$

for any choice of $P(a_{i},a_{j})$. The (linear) Perron-Frobenius theorem for
positive linear operator ensures the existence of a unique positive
eigenvector of $\Psi_{C}$ if the graph (4)

$(i,j)\leftrightarrow(k,\ell)\quad\text{if}\quad P(a_{i},a_{j})_{k,\ell}>0$
(4)

is connected.

##  Appendix D Proof of Proposition 3

Let $a_{k}$ and $a_{l}$ be two histograms from $A$, with supports associated
respectively to blocks $(i,i)$ and $(j,j)$. Then any transport plan $P$ in
$\Pi_{a_{k},a_{l}}$ is zero outside of the block $(i,j)$, so $\langle
C,\,P\rangle=c_{i,j}\times|\\!|P|\\!|_{1}=c_{i,j}$. Hence $\Phi_{A}(C)=C$.

##  Appendix E Proof of Proposition 6

The map $\Phi_{A}^{\infty}$ is

$\Phi_{A}^{\infty}(C)=(\langle-C(a_{i}-a_{j}),\,a_{i}-a_{j}\rangle=\langle
C,\,\mathcal{A}_{i,j}\rangle)_{i,j}\in\mathbb{R}^{m\times m}$
$\quad\text{where}\quad\mathcal{A}_{i,j}\triangleq-(a_{i}-a_{j})(a_{i}-a_{j})^{\top}\in\mathbb{R}^{n\times
n}$

We define the operator mapping correlation kernel to Euclidean distances

$\Delta(K)\triangleq-(K+K^{\top})+\operatorname{diag}(K)\mathds{1}_{n}^{\top}+\mathds{1}_{n}\operatorname{diag}(K)^{\top}.$

One has the convenient formula

$\Phi_{A}^{\infty}(C)=-\Delta(A^{\top}CA).$

Note that
$\ker(\Phi_{A}^{\infty})=\left\\{C\;;\;C=-C^{\top}\right\\}\cup\\{a\mathds{1}_{n}^{\top}+\mathds{1}_{n}b^{\top}\\}$,
$\operatorname{rank}(\Phi_{A}^{\infty})=n(n-1)/2$.

If $u\in\mathbb{R}^{n}$ is an eigenvector of $\tilde{A}^{\top}=JA^{\top}$ with
eigenvalue $\lambda\in\mathbb{R}$, then

$\displaystyle\Phi_{A}^{\infty}(\Delta(uu^{\top}))=-2\Phi_{A}^{\infty}(uu^{\top})\quad\text{(linearity
and }\ker\Phi_{A}^{\infty})$
$\displaystyle\quad=2\Delta(A^{\top}uu^{\top}A)\text{ (because
}\Delta(JKJ)=\Delta(K))$
$\displaystyle\quad=2\Delta(JA^{\top}uu^{\top}AJ)=2\lambda^{2}\Delta(uu^{\top})$

If $\lambda$ is complex, there is a pair $(v_{1},v_{2})$ associated to
eigenvalues $(\lambda,\bar{\lambda})$. We want to show that
$C\triangleq\Delta([v_{1}v_{2}][v_{1}v_{2}]^{\top}$ is an eigenvector. One has
similarely to the previous computation

$\displaystyle\Phi_{A}^{\infty}(C))$
$\displaystyle=-2\Phi_{A}^{\infty}([v_{1}v_{2}][v_{1}v_{2}]^{\top})$
$\displaystyle=-2\Delta(A^{\top}[v_{1}v_{2}][v_{1}v_{2}]^{\top}A)$
$\displaystyle=2\Delta(JA^{\top}[v_{1}v_{2}][v_{1}v_{2}]^{\top}AJ)$
$\displaystyle=2\Delta([\lambda v_{1}\bar{\lambda}v_{2}][\lambda
v_{1}\bar{\lambda}v_{2}]^{\top})$
$\displaystyle=2(\lambda^{2}+\bar{\lambda}^{2})\Delta([v_{1}v_{2}][v_{1}v_{2}]^{T})$