11institutetext: University of Toulouse, IRIT/INP-ENSEEIHT, F-31071 Toulouse,
France
11email: {Florentin.Coeurdoux<EMAIL_ADDRESS>
22institutetext: Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRIStAL, F-59000
Lille, France
22email<EMAIL_ADDRESS>
# Learning Optimal Transport Between two Empirical Distributions with
Normalizing Flows††thanks: This work was supported by the Artificial Natural
Intelligence Toulouse Institute (ANITI, ANR-19-PI3A-0004), the AI Sherlock
Chair (ANR-20-CHIA-0031-01), the ULNE national future investment programme
(ANR-16-IDEX-0004) and the Hauts-de-France Region.

Florentin Coeurdoux✉ 11 0000-0002-2100-1438 Nicolas Dobigeon 11
0000-0001-8127-350X Pierre Chainais 22 0000-0003-4377-7584

###### Abstract

Optimal transport (OT) provides effective tools for comparing and mapping
probability measures. We propose to leverage the flexibility of neural
networks to learn an approximate optimal transport map. More precisely, we
present a new and original method to address the problem of transporting a
finite set of samples associated with a first underlying unknown distribution
towards another finite set of samples drawn from another unknown distribution.
We show that a particular instance of invertible neural networks, namely the
normalizing flows, can be used to approximate the solution of this OT problem
between a pair of empirical distributions. To this aim, we propose to relax
the Monge formulation of OT by replacing the equality constraint on the push-
forward measure by the minimization of the corresponding Wasserstein distance.
The push-forward operator to be retrieved is then restricted to be a
normalizing flow which is trained by optimizing the resulting cost function.
This approach allows the transport map to be discretized as a composition of
functions. Each of these functions is associated to one sub-flow of the
network, whose output provides intermediate steps of the transport between the
original and target measures. This discretization yields also a set of
intermediate barycenters between the two measures of interest. Experiments
conducted on toy examples as well as a challenging task of unsupervised
translation demonstrate the interest of the proposed method. Finally, some
experiments show that the proposed approach leads to a good approximation of
the true OT.

###### Keywords:

Normalizing flows Optimal transport Generative Model.

##  1 Introduction

The optimal transport (OT) problem was initially formulated by the French
mathematician Gaspard Monge. In his seminal paper published in 1781 [18], he
raised the following question: how to move a pile of sand to a target location
with the least possible effort or cost? The objective was to find the best way
to minimize this cost by a transport plan, without having to list all the
possible matches between the starting and ending points. More recently, thanks
to recent advances related to computational issues [22], OT has founded
notable successes with respect to applications ranging from image processing
and computer vision [21] to machine learning [2] and domain adaptation [7].

Normalizing flows (NFs) have also attracted a lot of interest in the machine
learning community, motivated in particular by their ability to model high
dimensional data [19, 15]. These deep generative models are characterized by
an invertible operator that associates any input data distribution with a
target distribution that is usually chosen to be Gaussian. They have the great
advantage of leading to tractable distributions, which eases direct sampling
and density estimation. Applications of these generative models include image
generation with real-valued non-volume preserving transformations (RealNVP)
[10] or generative flows using an invertible 1x1 convolution (GLOW) [14].

Motivated by the similarities between the problem of OT and the training of
NF, this paper proposes a neural architecture and a corresponding training
strategy that permits to learn an approximate Monge map between any two
empirical distributions. The proposed framework is based on a relaxation of
the Monge formulation of OT. To adapt the training loss to the flow-based
structure of the network, this loss function is supplemented with a Sobolev
regularisation to promote minimal efforts achieved by each flow. Numerical
simulations show that this regularisation results in a smoother and more
efficient trajectory. Interestingly, the discretization inherent to the flow-
based structure of the network implicitly provides intermediate transports
and, at the same time, Wasserstein barycenters [1]. To the best of our
knowledge, this is the first time that NFs are considered to address OT and
Wasserstein barycenter computation, up to interesting dimensions.

Contributions. Our contributions are twofold: i) Section 2 recalls the Monge
formulation of OT and proposes a relaxation in the case of a transport between
two empirical distributions. ii) Section 3 presents the generic framework
based on NFs and describes a particular instance to solve the OT problem.
Section 4 presents some experimental results illustrating the performance of
the proposed method. Section 6 concludes this paper.

##  2 Relaxation of the optimal transport problem

Let $\mu$ and $\nu$ be two probability measures with finite second order
moments. More general measures, for example on $\mathcal{X}=\mathbb{R}^{d}$
(where $d\in\mathbb{N}^{*}$ is the dimension), can have a density
$\mathrm{d}\mu(x)=p_{X}(x)dx$ with respect to the Lebesgue measure, often
noted $p_{X}=\frac{\mathrm{d}\mu}{\mathrm{d}x}$, which means that

$\forall
h\in\mathcal{C}\left(\mathbb{R}^{d}\right),\quad\int_{\mathbb{R}^{d}}h(x)\mathrm{d}\mu(x)=\int_{\mathbb{R}^{d}}h(x)p_{X}(x)\mathrm{d}x$
(1)

where $\mathcal{C}(\cdot)$ is the class of continuous functions. In the
remainder of this paper, $\mathrm{d}\mu(x)$ and $p_{X}(x)\mathrm{d}x$ will be
used interchangeably.

###  2.1 Background on optimal transport

Let consider $\mathcal{X}$ and $\mathcal{Y}$ two separable metric spaces. Any
measurable application $T:\mathcal{X}\rightarrow\mathcal{Y}$ can be extended
to the so-called push-forward operator $T_{\sharp}$ which moves a probability
measure on $\mathcal{X}$ to a new probability measure on $\mathcal{Y}$. For
any measure $\mu$ on $\mathcal{X}$, one defines the image measure $\nu=$
$T_{\sharp}\mu$ on $\mathcal{Y}$ such that

$\forall
h\in\mathcal{C}(\mathcal{Y}),\quad\int_{\mathcal{Y}}h(y)\mathrm{d}\nu(y)=\int_{\mathcal{X}}h(T(x))\mathrm{d}\mu(x).$
(2)

Intuitively, the application $T:\mathcal{X}\rightarrow\mathcal{Y}$ can be
interpreted as a function moving a single point from one measurable space to
another [22]. The operator $T_{\sharp}$ pushes each elementary mass of a
measure $\mu$ on $\mathcal{X}$ by applying the function $T$ to obtain an
elementary mass in $\mathcal{Y}$. The problem of OT as formulated by Monge is
now stated in a general framework. For a given cost function
$c:\mathcal{X}\times\mathcal{Y}\rightarrow[0,+\infty]$, the measurable
application $T:\mathcal{X}\rightarrow\mathcal{Y}$ is called the OT map from a
measure $\mu$ to the image measure $\nu=T_{\\#}\mu$ if it reaches the infimum

$\inf_{T}\left\\{\int_{\mathcal{X}}c(x,T(x))\mathrm{d}\mu(x):T_{\sharp}\mu=\nu\right\\}.$
(3)

Alternatively the Kantorovitch formulation of OT results from a convex
relaxation of the Monge problem (3). By defining $\Pi$ as the set of all
probabilistic couplings with marginals $\mu$ and $\nu$, it yields the optimal
$\pi$ that reaches

$\min_{\pi\in\Pi}\int_{\mathcal{X}\times\mathcal{Y}}c\left(\mathbf{x},\mathbf{y}\right)d\pi\left(\mathbf{x},\mathbf{y}\right)$
(4)

Under this formulation, the optimal $\pi$, which is a joint probability
measure with marginals $\mu$ and $\nu$, can be interpreted as the optimal
transportation map. It allows the Wasserstein distance of order $p$ between
$\mu$ and $\nu$ to be defined as

$W_{p}\left(\mu,\nu\right)\stackrel{{\scriptstyle\text{ def
}}}{{=}}\inf_{\pi\in\Pi}\left\\{\left(\underset{\begin{subarray}{c}\mathbf{x}\sim\mu\\\
\mathbf{y}\sim\nu\end{subarray}}{\mathbb{E}}d\left(\mathbf{x},\mathbf{y}\right)^{p}\right)^{\frac{1}{p}}\right\\}$
(5)

where $d(\cdot,\cdot)$ is a distance defining the cost function
$c\left(\mathbf{x},\mathbf{y}\right)=d\left(\mathbf{x},\mathbf{y}\right)^{p}$.
The Wasserstein distance is also known as the Earth mover’s distance. It
defines a metric over the space of square integrable probability measures.

###  2.2 Proposed relaxation of OT

OT boils down to a variational problem, i.e., it requires the minimization of
an integral criterion in a class of admissible functions. Given two
probability measures $\mu$ and $\nu$, the existence and uniqueness of an
operator $T$ that belongs to the class of bijective, continuous and
differentiable functions such that $T_{\sharp}\mu=\nu$ is not guaranteed. The
difficulty lies in the class defining these admissible functions. Indeed, even
when $\mu$ and $\nu$ are regular densities on regular subsets of
$\mathbb{R}^{d}$, the search for a transport map such that $T_{\sharp}\mu=\nu$
makes the problem (3) difficult in a general case. To overcome the difficulty
of solving this equation on $T_{\sharp}$, we propose to reformulate the
Monge’s OT statement by relaxing the equality on the operator defining the
image measure.

More precisely, the equality between the image measure $T_{\sharp}\mu$ and the
target measure $\nu$ is replaced by the minimization of their statistical
distance $d(T_{\sharp}\mu,\nu)$. The choice of the distance $d(\cdot,\cdot)$
is crucial because it determines the quality of the approximation of the image
measure by the transport map $T$. In this work, we propose to choose
$d(\cdot,\cdot)$ as the Wasserstein distance ${W}_{p}(\cdot,\cdot)$. This
choice will be motivated by the fact that this distance can be easily
approximated numerically without explicit knowledge of the probability
distributions $\mu$ and $\nu$, in particular when they are empirically
described by samples only. The relaxation of the Monge problem (3) can then be
written as

$\inf_{T}\left\\{{W}_{p}(T_{\sharp}\mu,\nu)+\lambda\int_{\mathcal{X}}c(x,T(x))\mathrm{d}\mu(x)\right\\}$
(6)

where the cost function defined in (3) is interpreted here as a regularisation
term adjusted by the hyperparameter $\lambda$.

######  Remark 1.

The relaxed formulation (6) relies on the Wasserstein distance between the
target measure $\nu$ and the image measure $T_{\sharp}\mu$. This term should
not be confused with the Wasserstein distance ${W}_{p}(\mu,\nu)$ which is the
infimum reached by the solution of the Kantorovitch’s formulation of OT (4).

###  2.3 Discrete formulation

In a machine learning context, the underlying continuous measures are
conventionally approximated by empirical point measures thanks to available
data samples. Therefore, in this paper, we are interested in discrete measures
and the empirical formulation of the OT problem. Within this framework, we
will consider $\mu$ and $\nu$ two discrete measures described by the
respective samples $\mathbf{x}=\left\\{x_{n}\right\\}_{n=1}^{N}$ and
$\mathbf{y}=\left\\{y_{n}\right\\}_{n=1}^{N}$ such that
${\mu}=\frac{1}{N}\sum_{n=1}^{N}\delta_{x_{n}}$ and
${\nu}=\frac{1}{N}\sum_{n=1}^{N}\delta_{y_{n}}$. In the following, an
empirical version of the criterion (6) is proposed in the case of discrete
measures.

The formulation (6) requires the evaluation of a Wasserstein distance whose
computation is not trivial in its original form, especially in high dimension.
An alternative consists in considering its rewriting in the form of the
_sliced-Wasserstein_ (SW) distance. The idea underlying the SW distance is to
represent a distribution defined in high dimension thanks to a set of
projected one-dimensional distributions for which the computation of the
Wasserstein distance is closed-form. Let $p_{X}$ and $p_{Y}$ denote the
probability distributions of the random variables $X$ and $Y$. For any vector
on the unit sphere $u\in\mathbb{S}^{d-1}$, the projection operator
$S_{u}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is defined as
$S_{u}(x)\triangleq\langle u,x\rangle$. The SW distance of order
$p\in[1,\infty)$ between $p_{X}$ and $p_{Y}$ can be written [4]

$SW_{p}\left(p_{X},p_{Y}\right)=\left(\int_{\mathbb{S}^{d-1}}{W}_{p}\left(S_{u\sharp}p_{X},S_{u\sharp}p_{Y}\right)^{p}du\right)^{\frac{1}{p}}$
(7)

where the distance $W_{p}(\cdot,\cdot)$ defining the integrand is now one-
dimensional, leading to an explicit computation by inversion of the cumulative
distribution functions. In the case where the distributions $p_{X}$ and
$p_{Y}$ are represented by the respective samples $\mathbf{x}$ and
$\mathbf{y}$, a numerical Monte Carlo approximation of the SW distance is

$\widehat{SW}_{p}(\mathbf{x},\mathbf{y})=\frac{1}{J}\sum_{j=1}^{J}{W}_{p}\left(\frac{1}{N}\sum_{n=1}^{N}\delta_{S_{u_{j}}({x}_{n})},\frac{1}{N}\sum_{n=1}^{N}\delta_{S_{u_{j}}({y}_{n})}\right)$
(8)

where $u_{1},\ldots,u_{J}$ are drawn uniformly on the sphere
$\mathbb{S}^{d-1}$. The empirical form of the relaxation of the Monge problem
(6) is then written as

$\min_{T}\left\\{\widehat{SW}_{p}(T(\mathbf{x}),\mathbf{y})+\lambda\sum_{n=1}^{N}c\left(x_{n},T\left(x_{n}\right)\right)\right\\}$
(9)

where, with a slight abuse of notations,
$T(\mathbf{x})\triangleq\\{T(x_{n})\\}_{n=1}^{N}$.

##  3 Normalizing flows to approximate OT

This section proposes to solve the problem (9) by restricting the class of the
operator $T$ to a class of invertible deep networks referred to as
normalisation flows. The structure and the main properties of these networks
are detailed in paragraph 3.1. The strategy proposed to train these networks
to solve the problem (9) is then detailed in paragraph 3.2.

###  3.1 Normalizing flows

Normalization flows are a flexible class of deep generative networks that
intend to learn a change of variable between two probability distributions
$p_{X}$ and $p_{Y}$ through an invertible transformation
$T_{\boldsymbol{\Theta}}:X\mapsto Y=T_{\boldsymbol{\Theta}}(X)$ parametrized
by $\boldsymbol{\Theta}$. In general, the distribution $p_{X}$ is only known
through samples $\mathbf{x}=\left\\{x_{n}\right\\}_{n=1}^{N}$ and, for
tractability purpose, the distribution $p_{Y}$ is chosen as a centered normal
distribution with unit variance. The parameters $\boldsymbol{\Theta}$ defining
the operator $T_{\boldsymbol{\Theta}}$ are then adjusted by maximizing the
likelihood associated with the observations $\mathbf{x}$ according to the
change of variable formula

$p_{X}(x)=p_{Y}\left(T_{\boldsymbol{\theta}}(x)\right)\left|\operatorname{det}J_{T_{\boldsymbol{\theta}}^{-1}}\right|$
(10)

with $J_{T_{\boldsymbol{\Theta}}^{-1}}=\frac{\partial
T_{\boldsymbol{\theta}}^{-1}}{\partial x}$. NF networks obey a cell-like
structure, explicitly defining the operator $T_{\boldsymbol{\Theta}}(\cdot)$
as the composition of $M$ functions $T^{(m)}_{\boldsymbol{\theta}_{m}}$,
usually referred to as _flows_ , i.e.,

$T_{\boldsymbol{\Theta}}(\cdot)=T^{(M)}_{\boldsymbol{\theta}_{M}}\circ
T^{(M-1)}_{\boldsymbol{\theta}_{M-1}}\circ\ldots\circ
T^{(1)}_{\boldsymbol{\theta}_{1}}(\cdot)$ (11)

with
${\boldsymbol{\Theta}}=\left\\{\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{M}\right\\}$.
In the following, to lighten notations, each sub-function contributing to the
flow will be denoted by $T_{m}=T^{(m)}_{\boldsymbol{\theta}_{m}}$. In the
present work, these functions are chosen as coupling layers as implemented by
flows like RealNVP [10] and nonlinear independent component estimation (NICE)
[9]. These coupling layers ensure an invertible transformation and an explicit
expression of the Jacobian required in the change of variables (10). The input
and output of the $m$th layer are related as $\left(y_{\text{id
}},y_{\text{ch}}\right)=T_{m}\left(x_{\text{id }},x_{\text{ch}}\right)$ with

$\begin{cases}y_{\text{id }}=x_{\text{id }}\\\
y_{\text{ch}}=\left(x_{\text{ch}}+\mathsf{D}_{m}\left(x_{\text{id
}}\right)\right)\odot\exp\left(\mathsf{E}_{m}\left(x_{\text{id
}}\right)\right)\end{cases}$ (12)

where $x_{\text{id }}$ and $x_{\text{ch}}$ (resp. $y_{\text{id }}$ and
$y_{\text{ch}}$) are disjoint subsets of components of the input vector $x$
(resp. the output vector $y$). The splitting of the input $x$ into
$x_{\mathrm{id}}$ and $x_{\text{ch}}$ is achieved by a masking process such
that $x_{\text{ch}}=\textrm{mask}(x)$ is transformed into a function of the
unchanged part $x_{\text{id }}$. The scale function $\mathsf{E}_{m}(\cdot)$
and the offset function $\mathsf{D}_{m}(\cdot)$ are then described by neural
networks whose parameters $\boldsymbol{\theta}_{m}$ need to be adjusted during
the training. It is worth noting that imposing the flow-based architecture
detailed in (11) will lead to an explicit discretization scheme of the
transport map $T_{\boldsymbol{\Theta}}(\cdot)$ into a sequence of elementary
transport functions $T_{m}(\cdot)$. As it will be shown in Section 3.3, this
discretization has the great advantage of providing Wasserstein barycenters
associated with the two measures $\mu$ and $\nu$. Note that the proposed
method is not limited to NFs composed of coupling layers such as RealNVP [10],
NICE [9] or GLOW [14]. It can be generalized to other types of NFs, including
free-form Jacobian of reversible dynamics (FFJORD) [12] and masked
autoregressive flows (MAF) [20].

###  3.2 Loss function

As mentioned before, the objective of this work is to learn a bijective
operator relating any two distributions $p_{X}$ and $p_{Y}$ described by
samples $\mathbf{x}$ and $\mathbf{y}$. The search for this operator is
restricted to the class of invertible deep networks $T_{\boldsymbol{\Theta}}$
described in paragraph 3.1. The conventional strategy to train the network
would be to maximize the likelihood defined by (10). However this approach
cannot be implemented in the context of interest here since the base
distribution $p_{Y}$ is no longer explicitly given: it is only available
through the knowledge of the set of samples $\mathbf{y}$. As a consequence, to
adjust the weights of the network, the proposed alternative interprets the
underlying learning task as the search for a transport map. Then a first idea
would be to adjust these weights by directly solving the problem (9). However,
to take advantage of the flow-based architecture of the operator
$T_{\boldsymbol{\Theta}}(\cdot)$, it seems legitimate to equally distribute
the transport efforts provided by each flow. Thus, the regularization in (9)
will be instantiated for each elementary transformation $T_{m}(\cdot)$
associated to each flow of the network.

Moreover, when fitting deep learning-based models a major challenge arises
from the stochastic nature of the optimization procedure, which imposes to use
partial information (e.g., as mini-batches) to infer the whole structure of
the optimization landscape. On top of that, the cost function to be optimized
is not numerically constant since the approximation $\widehat{SW}$ of the SW
distance in (9) depends on the precise set of random vectors
$\left\\{{u}_{j}\right\\}_{j=1}^{J}$ drawn over the unit sphere. To alleviate
these optimization difficulties, we propose to further regularize the
objective function by penalizing the energy
$\left|{J_{T_{m}}(\cdot)}\right|^{2}$ of the Jacobians associated with the
transformations $T_{m}(\cdot)$, $m=1,\ldots,M$. These Sobolev-like penalties
promote regular operators $T_{m}(\cdot)$, promoting an overall operator
$T_{\boldsymbol{\theta}}(\cdot)$ regular itself [13]. In the context of
optimal transport, this regularization has already been studied in depth in
[17]. In that work, the author focused on the penalization of the Monge’s
formulation of OT by the $\ell_{2}$-norm of the Jacobian. It stated the
existence of an optimal transport map $T$ solving the minimization problem

$\inf_{T}\left\\{\int_{\mathcal{X}}\left(|T(x)-x|^{2}+\gamma|J_{T}|^{2}\right)T(x)\mathrm{d}x:T_{\\#}\mu=\nu\right\\}$
(13)

This formulation of OT imposes the transport map $T$ to be regular rather than
deducing its regularity from its optimal properties. Finally, the training of
the NF is carried out by minimizing the loss function

$\underbrace{\widehat{SW}_{p}(\mathbf{x},\mathbf{y})}_{\text{SW}}+\underbrace{\sum_{n=1}^{N}\sum_{m=1}^{M}\bigg{[}\lambda
c(T_{m-1}(x_{n}),T_{m}(x_{n}))+\gamma\left|{J_{T_{m}}}(x_{n})\right|^{2}\bigg{]}}_{\text{Reg}}$
(14)

with $T_{0}(x_{n})=x_{n}$. The proposed network, whose general architecture is
depicted in Fig. 1, will be referred to as SWOT-Flow in what follows.

Figure 1: Architecture of the proposed SWOT-Flow.

###  3.3 Intermediate transports and Wasserstein barycenters

As a consequence of the multiple-flow architecture (11) of the NF, the
transport map operated by the proposed SWOT-Flow is a composition of the $M$
individual flows $T_{m}(\cdot)$ ($m=1,\ldots,M$). Thus each flow implements an
elementary transport and the composition of the first $m$ flows defined as

$T_{[m]}(\cdot)\triangleq T_{m}\circ\ldots\circ T_{1}(\cdot)$ (15)

can be interpreted as an intermediate step of the transport map from the input
measure $\mu$ towards the target measure $\nu$, with $T_{[M]}(\cdot)\triangleq
T_{\boldsymbol{\Theta}}(\cdot)$. Interestingly, these intermediate transports
can be related to Wasserstein barycenters between $\mu$ and $\nu$ defined by
[1]

$\inf_{\beta}\left\\{\alpha
W_{p}(\mu,\beta)+(1-\alpha)W_{p}(\beta,\nu)\right\\}.$ (16)

Indeed, the next section dedicated to numerical experiments will empirically
show that $T_{[m]\sharp}\mu$ approaches the solution of the problem (16) for
the specific choice of the weight $\alpha=\frac{m}{M}$. In other words, the
image measures provided by each intermediate transport operated by SWOT-Flow,
i.e., as the outputs of each of the $M$ flows, can legitimately be interpreted
as Wasserstein barycenters.

##  4 Numerical experiments

This section assesses the versatility and the accuracy of SWOT-Flow through
two sets of numerical experiments. First, several toy experiments are
presented to provide some insights about key ingredients of the proposed
approach. Then the performance of SWOT-Flow is illustrated through the more
realistic and challenging task of unsupervised alignment of word embeddings in
natural language processing. The source code is publicly available on GitHub
111FlorentinCDX/SWOT-Flow.

###  4.1 Toy examples

In these experiments, the proposed framework SWOT-Flow is implemented and
tested with synthetic data. In all experiments, the input distributions are
described by the respective samples
$\mathbf{x}=\left\\{x_{n}\right\\}_{n=1}^{N}$ and
$\mathbf{y}=\left\\{y_{n}\right\\}_{n=1}^{N}$ such that
${\mu}=\frac{1}{N}\sum_{n=1}^{N}\delta_{x_{n}}$ and
${\nu}=\frac{1}{N}\sum_{n=1}^{N}\delta_{y_{n}}$ with $N=20000$. The cost
function $c(\cdot,\cdot)$ is chosen as the squared Euclidean distance, i.e.,
$c(x,y)=\|x-y\|_{2}^{2}$. However, it is worth noting that the proposed method
is not limited to this Euclidean distance and can handle other costs defined
on $\mathbb{R}^{d}$ or even on curved domains.

####  4.1.1 Implementation details.

The stochastic gradient descent used to solve (14) is implemented in Pytorch.
We use Adam optimizer with learning rate $10^{-4}$ and a batch size of $4096$
or $8192$ samples. The NF implementing $T_{\boldsymbol{\Theta}}(\cdot)$ is a
RealNVP [10] for the example of Fig. 2 and an ActNorm type architecture
network [14] for Fig. 3 and Fig. 4. It is composed of $M=4$ flows, each
composed of two four-layer neural networks corresponding to
$\mathsf{D}_{m}(\cdot)$ and $\mathsf{E}_{m}(\cdot)$ ($d\rightarrow
8\rightarrow 8\rightarrow d$) using hyperbolic tangent activation function.
During training, the number $J$ of slices drawn to approximate the SW distance
in (8) has been progressively increased, starting from $J=500$ to $J=2000$ by
step of $50$ slices. At each epoch, new slices are uniformly drawn over the
unit sphere and 100 epochs are carried out for each number of slices. The
training procedure consist in 1) defining the loss function as the sole SW
term in (14) from $J=500$ to $1500$ slices and then 2) incorporating the
regularization term denoted as Reg in (14) where hyperparameters $\lambda$ and
$\gamma$ are increased by a factor of $5\%$ every step of 100 slices.

Figure 2: Operator $T$ learnt by SWOT-Flow when the base distribution $p_{X}$
is a double-moon (top left) and the target distribution $p_{Y}$ is a circle
(bottom right).

####  4.1.2 Qualitative results.

As a first illustration of the flexibility of the proposed approach, Fig. 2
shows the results obtained after learning an operator $T$ that transports a
double moon-shaped distribution $p_{X}$ (top left) to a circle-shaped
distribution (bottom right). The empirical image measures $T_{\sharp}p_{X}$
(top right) and $T^{-1}_{\sharp}p_{Y}$ (bottom left) are obtained by applying
the estimated $T(\cdot)$ operator or its inverse $T^{-1}(\cdot)$. It is worth
noting that the difficulty inherent to this experiment lies in the respective
disjoint and non-disjoint supports of the two distributions. Despite the
regularity of the trained NF, a very good approximation of the OT is learnt,
even in presence of this topological change.

Figure 3: Elementary transports achieved by the proposed NF when trained
without (1st and 3rd panels) or with (2nd and 4th panels) the regularization
for protocols $\mathcal{P}_{1}$ (left panels) and $\mathcal{P}_{2}$ (right
panels).

Fig. 3 aims at illustrating the relevance of the Sobolev-like regularization
(i.e., the $\ell_{2}$-norm of the Jacobian) included into the loss function
(14) defined to train the NF. The first simulation protocol considers circle-
shaped distributions while the second case considers rectangle-shaped
distributions. In what follows, these two cases will be referred to as
$\mathcal{P}_{1}$ and $\mathcal{P}_{2}$, respectively. In this experiment, the
objective is to learn the transport map from an initial distribution $p_{X}$
(light blue) to a target distribution $p_{Y}$ (dark blue) which is translated
for $\mathcal{P}_{1}$ and both translated and stretched for $\mathcal{P}_{2}$.
The color gradient shows the outputs of the $M$ successive flows of the
network, i.e. the image measures $T_{[m]\sharp}p_{X}$ for $m=1,\ldots,M$. In
the absence of regularization (left), the successive elementary transports
clearly suffer from multiple unexpected deformations (superfluous translations
and dilations). In contrast, when the loss is complemented with the proposed
Sobolev-type penalty (right), the learnt operator $T$ is decomposed as a
sequence of much more regular elementary transports. The resulting transport
appears to be very close to optimal. In case $\mathcal{P}_{1}$, the expected
translation is recovered, as well as the combined translation and stretching
in case $\mathcal{P}_{2}$.

Table 1: Overall cost $\bar{C}$ and elementary costs $\bar{c}_{m}$ required by each flow $T_{m}(\cdot)$ of the NF trained with or without (w/o) regularization for protocols $\mathcal{P}_{1}$ (circle-shaped distributions) and $\mathcal{P}_{2}$ (rectangle-shaped distributions). | $\bar{c}_{1}$ | $\bar{c}_{2}$ | $\bar{c}_{3}$ | $\bar{c}_{4}$ | $\bar{C}$
---|---|---|---|---|---
$\mathcal{P}_{1}$ | w/o regularization | 150.13 | 110.94 | 108.41 | 151.65 | 521.12
with regularization | 90.20 | 90.70 | 90.71 | 90.22 | 361.22
$\mathcal{P}_{2}$ | w/o regularization | 154.99 | 98.67 | 52.49 | 101.21 | 407.38
with regularization | 88.77 | 89.42 | 89.43 | 89.38 | 357.0

To be more precise quantitatively, Table 1 compares some metrics obtained when
the NF has been trained using the regularization-free or regularized loss
function, as defined in (14). For the two aforementioned simulation protocols,
it reports the elementary costs

$\bar{c}_{m}=\frac{1}{N}\sum_{n=1}^{N}\left\|T_{m-1}(x_{n})-T_{m}(x_{n})\right\|_{2}^{2}$
(17)

spent by each of the M flows $T_{1}(\cdot),\ldots,T_{M}(\cdot)$ to achieve the
transport maps retrieved by SWOT-Flow. This table (last column) also reports
the overall cost $\bar{C}=\sum_{m=1}^{M}\bar{c}_{m}$. For the two simulation
protocols $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$, these results clearly show
cheaper transports when using the proposed regularization. For instance, for
the simulation protocol $\mathcal{P}_{1}$, the overall cost is $\bar{C}=360$
with the regularization, compared to $\bar{C}=520$ when it is omitted.
Moreover, when using the regularized loss function, this cost is distributed
homogeneously over the successive flows, with a variation of at most $\pm 1\%$
from one flow to another, against $\pm 20\%$ otherwise.

(a) Training set

(b) Squares

(c) Triangles

(d) Gaussians

(e) Training set

(f) Circles

(g) Squares

(h) Gaussians

Figure 4: Examples of transported data sets for protocols $\mathcal{P}_{1}$
(top) and $\mathcal{P}_{1}$ (bottom).

Fig. 4 aims at illustrating the capacity of generalization of the transport
map learnt by SWOT-Flow. In this experiment, SWOT-Flow has been trained
following the simulation protocols $\mathcal{P}_{1}$ (Fig. 4(a)) or
$\mathcal{P}_{2}$ (Fig. 4(e)). Once trained on the data set associated with
each protocol, the NFs are fed with differently shaped data and the elementary
transports are monitored as above. Fig 4(b)-4(d) and 4(g)-4(h) show the
results when using square-, triangle-, Gaussian-shaped data sets for both
protocols, respectively. As expected, all initial distributions are either
simply translated in case $\mathcal{P}_{1}$ or translated and stretched in
case $\mathcal{P}_{2}$. The intermediate distributions correspond to the
expected barycenters as well. Fig. 4 clearly demonstrates the generalization
capacity of the proposed approach.

####  4.1.3 Multivariate Gaussians with varying dimensions.

When the source and target distributions $\mu$ and $\nu$ of a transportation
problem are multivariate Gaussians, the Wasserstein barycenters defined by
(16) are also multivariate Gaussian distributions. In this case, an efficient
fixed-point algorithm can be used to estimate its mean vector $\mathbf{a}$ and
covariance matrix $\boldsymbol{\Sigma}$ [11]. This experiment capitalizes on
this finding to assess the ability of SWOT-Flow to approximate Wasserstein
barycenters, as stated in Section 3.3. To this end, the algorithm designed in
[11] is implemented to estimate the actual barycenter associated with two
prescribed multivariate Gaussian distributions for
$\alpha=1-\alpha=\frac{1}{2}$. This barycenter is compared to the image
measure $T_{[m]\sharp}\mu$ estimated by SWOT-Flow with $m=\frac{M}{2}$. More
precisely, the mean vector and the covariance matrix of the barycenter are
compared to their maximum likelihood estimates $\hat{\mathbf{a}}$ and
$\hat{\boldsymbol{\Sigma}}$ computed from the samples
$\\{T_{[m]}(x_{n})\\}_{n=1}^{N}$ transported by the first $m$ flows. The
resulting mean square errors (MSEs)

$\displaystyle\mathrm{MSE}(\mathbf{a})=\left\|{\mathbf{a}}-\hat{\mathbf{a}}\right\|_{2}^{2}\quad\text{and}\quad\mathrm{MSE}(\boldsymbol{\Sigma})=\|{\boldsymbol{\Sigma}}-\hat{\boldsymbol{\Sigma}}\|_{\mathrm{F}}^{2}$
(18)

are reported in Table 2 for varying dimensions ranging from $2$ to $8$. This
table also reports the MSEs reached by other state-of-the-art free-support
methods [8, 5, 16]. For the methods [8] and [5], $n=5000$ and $n=100$ support
points have been used, respectively, since these are the maximum numbers
allowed for the algorithms to terminate in reasonable computational times.
SWOT-Flow compares favorably to state-of-the-art methods since reported MSEs
in Table 2 appear to be most often the smallest. These observation may call
for a more general study, but remains noticeable since SWOT-Flow has not been
specifically designed to compute the Wasserstein barycenters, contrary to
alternate methods.

Table 2: Performance of the estimation of the median barycenters. Reported scores result from the average over 5 Monte Carlo runs. | [8] | [5] | [16] | SWOT-Flow
---|---|---|---|---
Dimension | 2 | $\mathrm{MSE}(\mathbf{a})$ | $9.99\cdot 10^{-5}$ | $3.14\cdot 10^{-4}$ | $1.17\cdot 10^{-4}$ | $\mathbf{8.09\cdot 10^{-5}}$
$\mathrm{MSE}(\boldsymbol{\Sigma})$ | $7.28\cdot 10^{-4}$ | $2.39\cdot 10^{-3}$ | $1.98\cdot 10^{-3}$ | $\mathbf{1.44\cdot 10^{-4}}$
4 | $\mathrm{MSE}(\mathbf{a})$ | $1.73\cdot 10^{-3}$ | $1.68\cdot 10^{-3}$ | $1.44\cdot 10^{-3}$ | $\mathbf{1.44\cdot 10^{-4}}$
$\mathrm{MSE}(\boldsymbol{\Sigma})$ | $1.35\cdot 10^{-2}$ | $2.50\cdot 10^{-2}$ | $1.22\cdot 10^{-2}$ | $\mathbf{3.61\cdot 10^{-4}}$
6 | $\mathrm{MSE}(\mathbf{a})$ | $2.04\cdot 10^{-3}$ | $\mathbf{2.58\cdot 10^{-3}}$ | $3.24\cdot 10^{-3}$ | $1.23\cdot 10^{-2}$
$\mathrm{MSE}(\boldsymbol{\Sigma})$ | $4.38\cdot 10^{-2}$ | $8.86\cdot 10^{-2}$ | $2.37\cdot 10^{-2}$ | $\mathbf{5.29\cdot 10^{-4}}$
8 | $\mathrm{MSE}(\mathbf{a})$ | $\mathbf{1.23\cdot 10^{-3}}$ | $1.48\cdot 10^{-3}$ | $3.14\cdot 10^{-3}$ | $1.29\cdot 10^{-2}$
$\mathrm{MSE}(\boldsymbol{\Sigma})$ | $8.31\cdot 10^{-2}$ | $1.64\cdot 10^{-1}$ | $4.23\cdot 10^{-2}$ | $\mathbf{2.22\cdot 10^{-3}}$

###  4.2 Unsupervised word translation

In a second set of experiments, the performance of SWOT-Flow has been assessed
on the task of unsupervised word translation. Given word embeddings trained on
two monolingual corpora, the goal is to infer a bilingual dictionary by
aligning the corresponding word vectors.

####  4.2.1 Experiment description.

This experiment considers the task of aligning two sets of points in high
dimension. More precisely, it aims at inferring a bilingual lexicon, without
supervision, by aligning word embeddings trained on monolingual data. FastText
[3] has been implemented to learn the word vectors used for representation. It
provides monolingual embeddings of dimension 300 trained on Wikipedia corpora.
Words are lower-cased, and those that appear less than 5 times are discarded
for training. As a post-processing step, only the first $50$k most frequent
words are selected in the reported experiments.

####  4.2.2 Architecture.

The proposed SWOT-Flow method has been implemented using a RealNVP
architecture. The scale function $\mathsf{E}_{m}(\cdot)$ and the offset
function $\mathsf{D}_{m}(\cdot)$ are multilayer neural networks with two
hidden layers of size 512 and hyperbolic tangent activation function. Adam has
been used as an optimizer with a learning rate of $1\cdot 10^{-3}$. The number
of slices involved in the Monte Carlo approximation of the SW distance in (8)
has been progressively increased from $J=500$ slices to $J=3000$ by steps of
50. For each number of slices, 100 epochs have been performed. The
hyperparameters $\lambda$ and $\gamma$ adjusting the weights of the composite
regularization have been increased by a factor of $5\%$ every steps of 500
slices.

Table 3: Comparison of accuracies obtained by SWOT-Flow and adv-net [6] for unsupervised word translation (’en’ is English, ’fr’ is French, ’de’ is German, ’ru’ is Russian). Method | en-es | es-en | en-fr | fr-en | en-de | de-en | en-ru | ru-en
---|---|---|---|---|---|---|---|---
SWOT-Flow | 20-NN | $37.4$ | $24.2$ | $46.6$ | $34.1$ | $44.4$ | $27.6$ | $14.4$ | $3.8$
10-NN | $\mathbf{33.5}$ | $\mathbf{22.5}$ | $\mathbf{42.5}$ | $32.5$ | $39.5$ | $26.8$ | $\mathbf{10.2}$ | $2.1$
adv-net [6] | 10-NN | $31.4$ | $21.2$ | $39.6$ | $\mathbf{35.1}$ | $\mathbf{40.1}$ | $\mathbf{27.1}$ | $7.1$ | $\mathbf{2.3}$

####  4.2.3 Main results.

To quantitatively measure the quality of SWOT-Flow, the problem of bilingual
lexicon induction is addressed, with the same setting as in [6]. The same
evaluation data sets and codes, as well as the same word vectors have been
used. Given an input word embedding ($n=1,\ldots,N$ with
$N_{\textrm{test}}=1000$) in a given language, the objective is to assess if
its counterpart $T(x_{n})$ transported by SWOT-Flow belongs to the close
neighborhood of the output word embedding $y_{n}$ in the target language. The
neighborhood $\mathcal{V}(y_{n})$ is defined as the set of $K$-nearest
neighbors computed in a cosine similarity sense with $K=10$ or $20$ in
dimension 300. The overall accuracy is computed as the percentage of correctly
transported input samples. Denoting by $\boldsymbol{1}_{\mathsf{A}}$ the
indicator function, i.e., $\boldsymbol{1}_{\mathsf{A}}=1$ if the assertion
$\mathsf{A}$ is true and $\boldsymbol{1}_{\mathsf{A}}=0$ otherwise,

$\mathrm{accuracy}=\frac{1}{N_{\textrm{test}}}\sum_{n=1}^{N_{\textrm{test}}}\boldsymbol{1}_{\left\\{T(x_{n})\in\mathcal{V}(y_{n})\right\\}}\times
100\ \text{(\%)}$ (19)

Table 3 reports the accuracy scores for several pairs of languages. Although
SWOT-Flow has not been specifically designed to perform word translation,
these results show that its overall performance is on par with the adversarial
network (adv-net) proposed specifically for this task in [6]. In particular,
SWOT-Flow seems to perform well for translation between languages with close
origins.

(a) english $\rightarrow$ french

(b) english $\rightarrow$ german

Figure 5: 2D PCA representation of the target word embedding space: the
targeted translated (in green) and the transported source (in blue) embedded
words.

Fig. 5 qualitatively illustrates this good performance by showing how close a
set of translated words $T(x_{n})$ are to their true translation $y_{n}$. This
representation is obtained by a classical projection on the 2 first PCA
components of the target embedded space. The translation of 5 specific words
from English to French or German fall in the close vicinity of their true
counterparts.

##  5 Discussion

####  5.0.1 Cycle consistency.

Cycle consistency, as proposed in CycleGAN [23], aims at learning meaningful
cross-domain mappings such that the data translated from the domain
$\mathcal{X}$ to the domain $\mathcal{Y}$ via
$T_{\mathcal{X}\rightarrow\mathcal{Y}}$ can be mapped back to the original
data points in $\mathcal{X}$ via $T_{\mathcal{Y}\rightarrow\mathcal{X}}$. That
is, $T_{\mathcal{Y}\rightarrow\mathcal{X}}\circ
T_{\mathcal{X}\rightarrow\mathcal{Y}}(x)\approx x$ for all $x\in\mathcal{X}$.
For CycleGan, and many other domain transfer models such as [2], this key
property should be enforced by including a cycle consistency term into the
loss function. Conversely, since NF-based generative models learn bijective
mappings, the proposed SWOT-Flow inherits the cycle consistency property by
construction.

####  5.0.2 Semi-discrete formulation.

The proposed SWOT-Flow framework has been explicitly derived to approximate OT
between two discrete empirical distributions. It can be instanciated to
perform semi-discrete OT, i.e., to handle the case where one of distribution
is not described by data points but rather given as an explicit continuous
probability measure. Instead of relaxing the Monge formulation (3) as in (6),
it would consist in replacing the SW distance with a log-likelihood term $\log
f_{\nu}(\cdot)$ associated with the target continuous measure. The loss
function in (14) would be replaced by

$-\sum_{n=1}^{N}\log
f_{\nu}(T_{\boldsymbol{\Theta}}({x_{n}}))+\sum_{n=1}^{N}\sum_{m=1}^{M}\bigg{[}\lambda
c(T_{m-1}(x_{n}),T_{m}(x_{n}))+\gamma\left|{J_{T_{m}}}(x_{n})\right|^{2}\bigg{]}$
(20)

where the log-likelihood term is evaluated at the data points
$\left\\{T_{\boldsymbol{\Theta}}({x_{n}})\right\\}_{n=1}^{N}$ transported by
the NF.

####  5.0.3 NF to approximate barycenters.

As discussed in Section 3.3 and experimentally illustrated in Section 4.1, the
flow-based architecture of the SWOT-Flow network leads to intermediate
transports, that can be related to Wasserstein barycenters. On the toy
Gaussian example considered in Section 4.1, SWOT-Flow provides good
approximation of the barycenters, although it has not been specifically
designed to perform this task. If one is interested in devising a NF
approximating these barycenters, the definition (16) would lead to the
optimization problem

$\inf_{T}\left\\{\sum_{m=1}^{M}\alpha_{m}W_{p}(\mu,T_{[m]\sharp}\mu)+(1-\alpha_{m})W_{p}(T_{[m]\sharp}\mu,\nu)\right\\}$
(21)

with $\alpha_{m}=\frac{m}{M}$. When handling empirical measures described by
samples, the subsequent discretization would require to replace both terms
with Monte Carlo approximations (8) of the SW distances. However, this would
lead to a computationally demanding training procedure.

##  6 Conclusion

We propose a new method to learn the optimal transport map between two
empirical distributions from sets of available samples. To this aim, we write
a relaxed and penalized formulation of the Monge problem. This formulation is
used to build a loss function that balances between the cost of the transport
and the proximity in Wasserstein distance between the transported base
distribution and the target one. The proposed approach relies on normalizing
flows, a family of invertible neural networks. Up to our knowledge, this is
the first method that is able to learn such a generalizable transport
operator. As a side benefit, the multiple flow architecture of the proposed
network interestingly yields intermediate transports and Wasserstein
barycenters. The proposed method is illustrated by numerical experiments on
toy examples as well as an unsupervised word translation task. Future work
will aim at extending these results to high dimensional applications.

## References

  * [1] Agueh, M., Carlier, G.: Barycenters in the Wasserstein space. SIAM J. Mathematical Analysis 43(2), 904–924 (2011)
  * [2] de Bézenac, E., Ayed, I., Gallinari, P.: Cyclegan through the lens of (dynamical) optimal transport. In: Joint Eur. Conf. Machine Learning and Knowledge Discovery in Databases (ECML-PKDD) (2021)
  * [3] Bojanowski, P., Grave, E., Joulin, A., Mikolov, T.: Enriching Word Vectors with Subword Information. Transactions of the Association for Computational Linguistics (2017)
  * [4] Bonneel, N., Rabin, J., Peyré, G., Pfister, H.: Sliced and radon Wasserstein barycenters of measures. J. Math. Imag. Vision 51(1), 22–45 (2015)
  * [5] Claici, S., Chien, E., Solomon, J.M.: Stochastic Wasserstein barycenters. In: Proc. Int. Conf. Machine Learning (ICML) (2018)
  * [6] Conneau, A., Lample, G., Ranzato, M., Denoyer, L., Jégou, H.: Word translation without parallel data. In: Proc. IEEE Int. Conf. Learn. Represent. (ICLR) (2018)
  * [7] Courty, N., Flamary, R., Tuia, D.: Domain adaptation with regularized optimal transport. In: Joint Eur. Conf. Machine Learning and Knowledge Discovery in Databases (ECML-PKDD) (2014)
  * [8] Cuturi, M., Doucet, A.: Fast computation of wasserstein barycenters. In: Xing, E.P., Jebara, T. (eds.) Proc. Int. Conf. Machine Learning (ICML) (2014)
  * [9] Dinh, L., Krueger, D., Bengio, Y.: NICE: non-linear independent components estimation. In: Bengio, Y., LeCun, Y. (eds.) Proc. IEEE Int. Conf. Learn. Represent. (ICLR) (2015)
  * [10] Dinh, L., Sohl-Dickstein, J., Bengio, S.: Density estimation using real NVP. In: Proc. IEEE Int. Conf. Learn. Represent. (ICLR) (2017)
  * [11] Álvarez Esteban, P.C., del Barrio, E., Cuesta-Albertos, J., Matrán, C.: A fixed-point approach to barycenters in wasserstein space. Journal of Mathematical Analysis and Applications (2016)
  * [12] Grathwohl, W., Chen, R.T.Q., Bettencourt, J., Sutskever, I., Duvenaud, D.: Ffjord: Free-form continuous dynamics for scalable reversible generative models. Proc. IEEE Int. Conf. Learn. Represent. (ICLR) (2019)
  * [13] Hoffman, J., Roberts, D.A., Yaida, S.: Robust learning with jacobian regularization. arXiv (2020)
  * [14] Kingma, D.P., Dhariwal, P.: Glow: Generative flow with invertible 1x1 convolutions. In: Adv. in Neural Information Process. Systems (NeurIPS) (2018)
  * [15] Kobyzev, I., Prince, S.J., Brubaker, M.A.: Normalizing flows: An introduction and review of current methods. IEEE Trans. Patt. Anal. Mach. Intell. 43(11), 3964–3979 (2020)
  * [16] Korotin, A., Li, L., Solomon, J., Burnaev, E.: Continuous wasserstein-2 barycenter estimation without minimax optimization. In: Proc. IEEE Int. Conf. Learn. Represent. (ICLR) (2021)
  * [17] Louet, J.: Problèmes de transport optimal avec pénalisation en gradient. Ph.D. thesis, Université Paris-Sud, France (2014)
  * [18] Monge, G.: Mémoire sur la théorie des déblais et des remblais. Histoire de l’Académie Royale des Sciences de Paris (1781)
  * [19] Papamakarios, G., Nalisnick, E., Rezende, D.J., Mohamed, S., Lakshminarayanan, B.: Normalizing flows for probabilistic modeling and inference. J. Mach. Learning Research 22(57), 1–64 (2021)
  * [20] Papamakarios, G., Pavlakou, T., Murray, I.: Masked autoregressive flow for density estimation. In: Adv. in Neural Information Process. Systems (NeurIPS) (2017)
  * [21] Paulin, L., Bonneel, N., Coeurjolly, D., Iehl, J.C., Webanck, A., Desbrun, M., Ostromoukhov, V.: Sliced optimal transport sampling. ACM Trans. Graphics (Proc. SIGGRAPH) (2020)
  * [22] Peyré, G., Cuturi, M.: Computational optimal transport: with applications to data science. Foundations and Trends® in Machine Learning 11(5-6), 355–607 (2019)
  * [23] Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired image-to-image translation using cycle-consistent adversarial networks. In: Proc. IEEE Int. Conf. Computer Vision (ICCV) (2017)


Persistent current and low-field magnetic susceptibility in one-dimensional
mesoscopic rings: Effect of long-range hopping

Santanu K. Maiti1,2,∗

1Theoretical Condensed Matter Physics Division, Saha Institute of Nuclear
Physics,
1/AF, Bidhannagar, Kolkata-700 064, India
2Department of Physics, Narasinha Dutt College, 129, Belilious Road,
Howrah-711 101, India

Abstract

Persistent current and low-field magnetic susceptibility in single-channel
normal metal rings threaded by a magnetic flux $\phi$ are investigated within
the tight-binding framework considering long-range hopping of electrons in the
shortest path. The higher order hopping integrals try to reduce the effect of
disorder by delocalizing the energy eigenstates, and accordingly, current
amplitude in disordered rings becomes comparable to that of an ordered ring.
Our study of low-field magnetic susceptibility predicts that the sign of
persistent currents can be mentioned precisely in mesoscopic rings with fixed
number of electrons, even in the presence of impurity in the rings. For
perfect rings, low-field current shows only the diamagnetic sign irrespective
of the total number of electrons $N_{e}$. On the other hand, in disordered
rings it exhibits the diamagnetic and paramagnetic natures for the rings with
odd and even $N_{e}$ respectively.

PACS No.: 73.23.Ra; 73.23.-b; 71.23.-k; 73.20.Jc; 75.20.-g

Keywords: Persistent current; Magnetic susceptibility; Long-range hopping;
Finite temperature; Disorder.

∗Corresponding Author: Santanu K. Maiti

Electronic mail<EMAIL_ADDRESS>
##  1 Introduction

The phenomenon of persistent current in mesoscopic normal metal rings has
generated a lot of excitement as well as controversy over the past years. In a
pioneering work, Büttiker, Imry and Landauer [1] predicted that, even in the
presence of disorder, an isolated one-dimensional metallic ring threaded by a
magnetic flux $\phi$ can support an equilibrium persistent current with
periodicity $\phi_{0}=ch/e$, the elementary flux quantum. Later, the existence
of persistent current was further confirmed by several experiments [2, 3, 4,
5, 6, 7, 8]. However, these experiments yield many results those are not well-
understood theoretically even today [9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19, 20, 21, 22, 23, 24, 25]. The results of the single loop experiments are
significantly different from those for the ensemble of isolated loops.
Persistent currents with expected $\phi_{0}$ periodicity have been observed in
isolated single Au rings [2] and in a GaAs-AlGaAs ring [3]. Levy et al. [4]
found oscillations with period $\phi_{0}/2$ rather than $\phi_{0}$ in an
ensemble of $10^{7}$ independent Cu rings. Similar $\phi_{0}/2$ oscillations
were also reported for an ensemble of disconnected $10^{5}$ Ag rings [5] as
well as for an array of $10^{5}$ isolated GaAs-AlGaAs rings [6]. In an
experiment, Jariwala et al. [7] obtained both $\phi_{0}$ and $\phi_{0}/2$
periodic persistent currents for an array of thirty diffusive mesoscopic Au
rings. Except for the case of nearly ballistic GaAs-AlGaAs ring [3], all the
measured currents are in general one or two orders of magnitude larger than
those expected from theory. The diamagnetic response of the measured
$\phi_{0}/2$ oscillations of ensemble-averaged persistent currents near zero
magnetic field also contrasts with most predictions [15, 16].

Free electron theory predicts that, at zero temperature, an ordered one-
dimensional metallic ring threaded by a magnetic flux $\phi$ supports
persistent current with maximum amplitude $I_{0}=ev_{F}/L$, where $v_{F}$ is
the Fermi velocity of an electron and $L$ is the circumference of the ring.
Metals are intrinsically disordered which tends to decrease persistent
current, and calculations show that the magnitude of the currents reduces to
$I_{0}l/L$, where $l$ is the elastic mean free path of the electrons. This
expression remains valid even if one takes into account the finite width of
the ring by adding contributions from the transverse channels, since disorder
leads to a compensation between the channels [10, 11]. However, measurements
on single isolated mesoscopic rings [2, 3] detected $\phi_{0}$-periodic
persistent currents with amplitudes of the order of $I_{0}\sim ev_{F}/L$,
(close to the value for an ordered ring). Though theory seems to agree with
experiment [3] only when disorder is weak, the amplitudes of the currents in
single-isolated-diffusive gold rings [2] are two orders of magnitude larger
than the theoretical estimates. This discrepancy initiated intense theoretical
activity, and it is generally believed that the electron-electron correlation
plays an important role in the disordered rings [17, 18, 19], though the
physical origin behind this enhancement of persistent current is still
unclear.

In this article we investigate a detailed study of persistent current and low-
field magnetic susceptibility in single channel rings in the tight-binding
framework considering long-range hopping of electrons in the shortest path.
Our calculations show that in a disordered ring with higher order hopping
integrals, current amplitude is comparable to that of an ordered ring. This is
due to the fact that higher order hopping integrals try to delocalize the
energy eigenstates and thus compensate the effect of disorder. In the rest
part of this article, we describe the dependences of the sign of low-field
currents as a function of the total number of electrons $N_{e}$, and also
discuss the effect of temperature on these low-field currents.

The plan of the paper is as follow. Section $2$ relates the behavior of
persistent current in the presence of long-range hopping integrals and clearly
describes how current amplitude in disordered rings becomes comparable to that
of an ordered ring. In Section $3$, we investigate the behavior of low-field
magnetic susceptibility at absolute zero temperature both for ordered and
disordered rings as a function of $N_{e}$. Section $4$ focuses the effect of
temperature on the low-field currents and determines the critical value of
magnetic flux $\phi_{c}(T)$ where current changes its sign from the
paramagnetic to diamagnetic phase for the rings with even number of electrons.
Finally, the conclusions of our study can be found in Section $5$.

##  2 Persistent current

We describe a $N$-site ring (see Fig. 1) enclosing a magnetic flux $\phi$ (in
units of the elementary flux quantum $\phi_{0}=ch/e$) by the following tight-
binding Hamiltonian in the Wannier basis,

$H=\sum_{i}\epsilon_{i}c_{i}^{\dagger}c_{i}+\sum_{i\neq
j}v_{ij}\left[e^{-i\theta}c_{i}^{\dagger}c_{j}+h.c.\right]$ (1)

where $\epsilon_{i}$’s are the site energies, $v_{ij}$’s are the hopping
integrals between the sites $i$ and $j$, and $\theta={2\pi\phi}(|i-j|)/N$. The
long-range hopping (LRH) integrals are taken as,

$v_{ij}=\frac{v}{\left|\frac{N}{\pi}\sin\left[\frac{\pi}{N}(i-j)\right]\right|^{\alpha}}$
(2)

where, $v$ being a constant representing the nearest-neighbor hopping (NNH)
integral. In the present work, electron-electron interaction is not included,
and therefore, we do not consider the spin of the electrons since it will not
change any qualitative behavior of persistent currents. Throughout our study
we set $v=-1$, and use the units where $c=e=h=1$.

An electron in an eigenstate with energy $E_{n}$ carries a persistent

Figure 1: One-dimensional mesoscopic normal metal ring threaded by a magnetic
flux $\phi$. A persistent current $I$ is established in the ring.

current $I_{n}(\phi)=-\partial E_{n}/\partial\phi$, and at zero temperature
total persistent current is given by $I(\phi)=\sum_{n}I_{n}(\phi)$, where the
summation is over all the states below the Fermi level.

For an ordered ring, setting $\epsilon_{i}=0$ for all $i$, the energy of the
$n$-th eigenstate can be expressed as,

$E_{n}(\phi)=\sum_{m}\frac{2v}{m^{\alpha}}\cos\left[\frac{2\pi
m}{N}(n+\phi)\right]$ (3)

where $m$ is an integer and it runs from $1$ to $N/2$ for the rings with even
$N$, while it goes from $1$ to $(N-1)/2$ for those rings described by odd $N$.
Now the persistent current carried by this $n$-th eigenstate becomes,

$I_{n}(\phi)=\left(\frac{4\pi
v}{N}\right)\sum_{m}m^{1-\alpha}\sin\left[\frac{2\pi m}{N}(n+\phi)\right]$ (4)

For very large value of $\alpha$, Eqs. (3) and (4) essentially reduce to the
expressions for the energy spectrum and persistent current of an ordered ring
described by the nearest-neighbor tight-binding Hamiltonian. As we decrease
$\alpha$, contributions from the higher order hopping integrals become
appreciable which modify the energy spectrum and persistent current, and we
will see that the modifications are quite significant in the presence of
disorder. Figure 2 shows the variation of persistent current as a function of
magnetic flux $\phi$ for some perfect rings ($N=120$). The solid and dotted
curves represent the results for the rings described by LRH ($\alpha=1.3$) and
NNH integrals respectively, where the curves plotted in (a) give the variation
of the currents for the rings with odd $N_{e}$

Figure 2: Current-flux characteristics for some perfect rings with size
$N=120$. The solid and dotted curves are respectively for the rings with LRH
($\alpha=1.3$) and NNH integrals, where (a) $N_{e}=35$ (odd) and (b)
$N_{e}=40$ (even).

($N_{e}=35$) and the same are plotted for the rings with even $N_{e}$
($N_{e}$=40) in (b). The results predict that the current amplitude increases
in the presence of LRH integrals, compared to the NNH integral. In this
context it has been examined that, in the presence of LRH integrals, the
amplitude of the current initially increases (not shown here in the figure) as
we increase the ring size, but eventually it falls when the ring becomes
larger. This is due to the fact that as we increase the number of sites, the
Hamiltonian Eq. (1) includes additional higher order hopping integrals which
cause an increase in the net velocity of the electrons, but after certain ring
size the increment in velocity drops to zero since the additional hopping
integrals are then between far enough sites giving negligible contributions.

Now we address the problem of persistent current in the presence of disorder.
In order to introduce the disorder in the ring, we choose the site energies
$\epsilon_{i}$’s randomly from a “Box” distribution function of width $W$,
which reveal that the ring is subjected to the diagonal disorder. As
representative examples of persistent current in disordered rings, we plot

Figure 3: Current-flux characteristics for some typical disordered rings
($W=1$) with size $N=120$. The solid and dotted curves correspond to the rings
with LRH ($\alpha=1.3$) and NNH integrals respectively, where (a) $N_{e}=35$
(odd) and (b) $N_{e}=40$ (even).

the results in Fig. 3 for some $120$-site rings taking $W=1$. All the curves
shown in Fig. 3 are performed for the distinct disordered configurations of
the rings and no averaging is done here since in the averaging process several
mesoscopic phenomena disappear. The solid and dotted curves correspond to the
rings with LRH and NNH integrals respectively, and our results predict that
current amplitude gets an order of magnitude enhancement in the rings
described by the LRH integrals compared to those rings described by the NNH
integral. Figure 3(a) shows the variation of the persistent current for the
rings with odd $N_{e}$ ($N_{e}=35$) and Fig. 3(b) gives the variation of the
currents for the rings with even $N_{e}$ ($N_{e}=40$). It is apparent from
Figs. 2 and 3 that the current amplitudes in disordered rings with LRH
integrals are of the same order of magnitude as observed in ordered rings. We
have also seen that the decrease in amplitude of the current is quite small
even if we increase the strength of disorder. In the NNH models current
amplitudes are suppressed due to the localization of the energy eigenstates
[26]. On the other hand, the present tight-binding model with LRH integrals
supports extended electronic eigenstates even in the presence of disorder and
for this reason persistent currents are not reduced by the impurities.

##  3 Magnetic susceptibility at $T=0$ K

The sign of persistent currents can be determined exactly by calculating
magnetic susceptibility, and here we investigate the properties of low-field

Figure 4: $\chi(\phi)$ versus $N_{e}$ curves for some perfect rings with
$N=150$. The solid and dotted curves represent the rings with LRH
($\alpha=1.4$) and NNH integrals respectively.

currents for the rings with fixed number of electrons $N_{e}$. The general
expression of magnetic susceptibility is expressed in the form,

$\chi(\phi)=\frac{N^{3}}{16\pi^{2}}\left(\frac{\partial
I(\phi)}{\partial\phi}\right)$ (5)

Thus by calculating the sign of $\chi(\phi)$ we can predict whether the
current is diamagnetic or paramagnetic.

Let us first discuss the properties of low-field currents in perfect rings at
absolute zero temperature ($T=0$ K). In Fig. 4, we plot $\chi(\phi)$ as a
function of $N_{e}$ in the limit $\phi\rightarrow 0$ for some perfect rings
with $N=150$. The solid and dotted curves represent the variation of $\chi$
for the rings described by LRH and NNH integrals respectively. The results
show that, for perfect rings, low-field currents exhibit only the diamagnetic
sign irrespective of the total number of electrons $N_{e}$ in the ring. This
variation can be clearly understood if we consider the slope of the curves as
plotted in Figs. 2(a) and (b)

Figure 5: $\chi(\phi)$ versus $N_{e}$ curves of some disordered rings ($W=1$)
with size $N=150$. The solid and dotted lines are respectively for the rings
with even and odd $N_{e}$, where (a) NNH and (b) LRH ($\alpha=1.4$) integrals.

near $\phi=0$. Thus, both for the perfect rings with odd and even $N_{e}$,
current has only negative slope which predicts the diamagnetic persistent
current.

The effect of impurities on the sign of low-field currents is quite
interesting. As representative examples, in Fig. 5 we display the variation of
$\chi$ as a function of $N_{e}$, where (a) and (b) represent the rings with
NNH and LNH integrals respectively. The solid and dotted lines correspond to
the results for the rings containing even and odd number of electrons
respectively. These results emphasize that, in the disordered rings the low-
field currents exhibit the diamagnetic sign for odd $N_{e}$, while we get the
paramagnetic response for the rings with even $N_{e}$. The diamagnetic and the
paramagnetic natures of the low-field currents in the presence of impurity in
the rings can be understood easily if we take the slope of the curves as given
in Figs. 3(a) and (b). Such an effect of disorder on the low-field currents is
true for any disordered configuration. Accordingly, in the presence of
impurity one can easily predict the sign of the low-field currents both for
the rings with odd and even $N_{e}$, irrespective of the specific realization
of disordered configuration of the rings.

##  4 Magnetic susceptibility at finite temperature

This section focuses the effect of temperature on the low-field currents.

Figure 6: Variation of $\phi_{c}(T)$ with $N_{e}$ (even $N_{e}$ only) for some
disordered rings ($W=1$) taking $N=60$, where (a) rings with NNH integrals and
(b) rings with LRH ($\alpha=1.6$) integrals. The upper and lower curves both
in (a) and (b) are respectively for the rings with $T/T^{\star}=1.0$ and
$T/T^{\star}=0.5$.

As temperature increases the probability that electrons occupy higher energy
levels, those may carry larger currents, increases. But if we increase the
temperature in such a way that crosses the energy gap between two successive
energy levels, which carry currents in opposite directions, then mutual
cancellations of the positive and negative currents decrease the net current
amplitude. Therefore, it is necessary to specify a characteristic temperature
$T^{\star}$, which is determined by the energy level spacing $\Delta$. At
finite temperature, thermal excitations, such as phonons, are present which
interact with the electrons inelastically and thus randomizes the phase of the
electronic wave functions. These interactions try to destroy the phase
coherence of the electrons which can remove the quantum effects. Hence, it is
necessary to do the calculations at sufficiently low temperatures such that
the phase coherence length of an electron exceeds the circumference $L$ of the
ring.

In our calculations of low-field magnetic susceptibility at absolute zero
temperature ($T=0$ K), we see that the current has only diamagnetic sign for
perfect rings irrespective of the total number of electrons, while in the
presence of impurity, it exhibits respectively the diamagnetic and
paramagnetic sign for the rings with odd and even $N_{e}$. It is well known
that at any finite temperature ($T\neq 0$ K), low-field current has a
paramagnetic sign for the rings with even $N_{e}$ and in this section we
determine that critical value of magnetic flux, $\phi_{c}(T)$, where the low-
field current changes its sign from the paramagnetic to diamagnetic nature. In
Fig. 6 we display the variation of $\phi_{c}(T)$ as a function of $N_{e}$
(here $N_{e}$ is even only) for $60$-site disordered rings ($W=1$) at two
different temperatures. The upper and lower curves in Figs. 6(a) and (b) are
respectively for the rings with $T/T^{\star}=1.0$ and $0.5$. Figure 6(a) shows
the results for the rings with NNH integral, while the same are plotted for
the rings with LRH ($\alpha=1.6$) integrals in Fig. 6(b). From these results
we can emphasize that, as the temperature increases the critical value of
magnetic flux $\phi_{c}(T)$, where the low-field current changes its sign from
the paramagnetic phase to the diamagnetic one, increases.

##  5 Concluding remarks

In conclusion, we have investigated the behavior of persistent current in
single-isolated mesoscopic rings subjected to both NNH and LRH integrals
within the tight-binding framework. Our exact numerical calculations have
shown that the current amplitude in disordered rings are comparable to that of
ordered rings if we consider the model with LRH integrals instead of usual NNH
integral models. This is due to the fact that higher order hopping integrals
try to delocalize the energy eigenstates and thus prevents the reduction of
current due to disorder in the rings. Later, we have studied the low-field
magnetic response at $T=0$ K both for the perfect and disordered rings and our
results have predicted that the sign of the currents can be mentioned
precisely even in the presence of impurity in the rings. At the end, we have
calculated the magnetic response at finite temperatures ($T\neq 0$ K) and
estimated the critical value of magnetic flux $\phi_{c}(T)$ where the low-
field current changes its sign from the paramagnetic to the diamagnetic
nature.

## References

  * [1] M. Büttiker, Y. Imry, and R. Landauer, Phys. Lett. 96A, 365 (1983).
  * [2] V. Chandrasekhar, R. A. Webb, M. J. Brady, M. B. Ketchen, W. J. Gallagher, and A. Kleinsasser, Phys. Rev. Lett. 67, 3578 (1991).
  * [3] D. Mailly, C. Chapelier, and A. Benoit, Phys. Rev. Lett. 70, 2020 (1993).
  * [4] L. P. Levy, G. Dolan, J. Dunsmuir, and H. Bouchiat, Phys. Rev. Lett. 64, 2074 (1990).
  * [5] R. Deblock, R. Bel, B. Reulet, H. Bouchiat, and D. Mailly, Phys. Rev. Lett. 89, 206803 (2002).
  * [6] B. Reulet, M. Ramin, H. Bouchiat, and D. Mailly, Phys. Rev. Lett. 75, 124 (1995).
  * [7] E. M. Q. Jariwala, P. Mohanty, M. B. Ketchen, and R. A. Webb, Phys. Rev. Lett. 86, 1594 (2001).
  * [8] W. Rabaud, L. Saminadayar, D. Mailly, K. Hasselbach, A. Benoit and B. Etienne, Phys. Rev. Lett. 86, 3124 (2001).
  * [9] H. F. Cheung, Y. Gefen, E. K. Riedel, and W. H. Shih, Phys. Rev. B 37, 6050 (1988).
  * [10] H. F. Cheung and E. K. Riedel, Phys. Rev. Lett. 62, 587 (1989).
  * [11] G. Montambaux, H. Bouchiat, D. Sigeti, and R. Friesner, Phys. Rev. B 42, 7647 (1990).
  * [12] H. Bouchiat and G. Montambaux, J. Phys. (Paris) 50, 2695 (1989).
  * [13] B. L. Altshuler, Y. Gefen, and Y. Imry, Phys. Rev. Lett. 66, 88 (1991).
  * [14] F. von Oppen and E. K. Riedel, Phys. Rev. Lett. 66, 84 (1991).
  * [15] A. Schmid, Phys. Rev. Lett. 66, 80 (1991).
  * [16] V. Ambegaokar and U. Eckern, Phys. Rev. Lett. 65, 381 (1990).
  * [17] M. Abraham and R. Berkovits, Phys. Rev. Lett. 70, 1509 (1993).
  * [18] G. Bouzerar, D. Poilblanc, and G. Montambaux, Phys. Rev. B 49, 8258 (1994).
  * [19] T. Giamarchi and B. S. Shastry, Phys. Rev. B 51, 10915 (1995).
  * [20] X. Waintal, G. Fleury, K. Kazymyrenko, M. Houzet, P. Schmitteckert and D. Weinmann, Phys. Rev. Lett. 101, 106804 (2008).
  * [21] K. Yakubo, Y. Avishai and D. Cohen, Phys. Rev. B 67, 125319 (2003).
  * [22] E. H. M. Ferreira, M. C. Nemes, M. D. Sampaio and H. A. Weidenmüller, Phys. Lett. A 333, 146 (2004).
  * [23] P. A. Orellana, M. L. Ladron de Guevara, M. Pacheco and A. Latge, Phys. Rev. B 68, 195321 (2003).
  * [24] I. O. Kulik, Physica B 284, 1880 (2000).
  * [25] M. V. Moskalets, Physica B 291, 350 (2000).
  * [26] P. A. Lee and T. V. Ramakrishnan, Rev. Mod. Phys. 57, 287 (1985).