# Adaptive Mesh Methods on Compact Manifolds via Optimal Transport and Optimal
Information Transport

Axel G. R. Turnquist  Department of Mathematical Sciences, New Jersey
Institute of Technology, University Heights, Newark, NJ 07102<EMAIL_ADDRESS>
###### Abstract.

Moving mesh methods are designed to redistribute a mesh in a regular way. This
applied problem can be considered to overlap with the problem of finding a
diffeomorphic mapping between density measures. In applications, an off-the-
shelf grid needs to be restructured to have higher grid density in some
regions than others. This should be done in a way that avoids tangling, hence,
the attractiveness of diffeomorphic mapping techniques. For exact
diffeomorphic mapping on the sphere a major tool used is Optimal Transport,
which allows for diffeomorphic mapping between even non-continuous source and
target densities. However, recently Optimal Information Transport was
rigorously developed allowing for exact and inexact diffeomorphic mapping and
the solving of a simpler partial differential equation.

In this manuscript, we solve adaptive mesh problems using Optimal Transport
and Optimal Information Transport on the sphere and introduce how to
generalize these computations to more general manifolds. We choose to perform
this comparison with provably convergent solvers, which is generally
challenging for either problem due to the lack of boundary conditions and lack
of comparison principle in the partial differential equation formulation.

The author was partially supported by an NSF GRFP 1849508 and by NSF DMS
1751996

Diffeomorphic density matching has found widespread applications in medical
image registration see, for example, [12, 13, 7, 27], random sampling from a
probability measure see, for example, [2, 24] and adaptive moving mesh
methods, for example, in meteorology see the line of work in [3, 31, 4], among
other applications. These applications are necessarily more challenging from a
theoretical perspective as well as a computational perspective when the
geometry is non-Euclidean. In this manuscript, we focus on the non-Euclidean
case. Diffeomorphic density matching techniques largely fall into three
categories: (1) classical methods, (2) Optimal Transport methods, and (3)
Optimal Information Transport methods.

Classical methods in imaging sciences consist of positive scalar image
functions composed from the right by transformations [33]. That is, given a
mapping $T$, the source density by a density function $f_{0}$ and the target
by $f_{1}$, classical techniques compute $f_{0}\circ T=f_{1}$. Typically,
classical methods consists of feature detection, feature matching, transform
model estimation, and image resampling and transformation, as explained in
[32]. Features can consist of Fourier components, eigenvalues, et cetera. Non-
classical methods, like Optimal Transport and Optimal Information Transport,
allow the transformation to act as a pullback on the density, that is
$\left|DT\right|f_{0}\circ T=f_{1}$. This generalization has particular
benefit in that it allows for proving the existence of a diffeomorphic mapping
over a much wider range of densities $f_{0}$ and $f_{1}$.

Optimal Transport density matching in non-Euclidean geometries is a field of
study in which much progress on the theoretical side has been made recently,
see the following for examples of this large line of work [22, 18, 19, 9, 29,
20]. Recently, the authors in [14, 15] introduced a numerical convergence
framework and convergent numerical scheme for the PDE formulation of the
optimal transport problem on the sphere which included the traditional squared
geodesic cost $c(x,y)=\frac{1}{2}d_{\mathbb{S}^{2}}(x,y)^{2}$. The convergence
framework addressed situations where the solutions $u$ of the PDE formulation
are known to be a priori smooth (that is, at least $C^{2}$) as well as
nonsmooth (but still $C^{1}$). This convergence framework required that the
numerical schemes were consistent, stable in the average value, and monotone
and achieved compactness via a Lipschitz constraint imposed on the PDE level.

Optimal Information Transport was developed in the papers [1, 23], which treat
both exact and inexact diffeomorphic density matching based on a principal
fiber bundle structure and the particularly nice way that horizontal geodesics
in an information metric are pulled down to geodesics with respect to the
Fisher-Rao metric. This allows for the exact diffeomorphic map between two
smooth densities to be computed via solving the Poisson equation. Limitations
to this approach are addressed as well, such as its apparent inappropriateness
for image interpolation problems, by introducing computations with a non-
compatible metric (which allows for the source density to vanish) [1] in
Euclidean space. Optimal Information Transport methods have subsequently found
applications in image registration for medical applications, for example in
[7, 27] and also random sampling [2], also in subsets of Euclidean space.

In this manuscript, we establish a class of provably convergent methods (with
convergence rate guarantees) for computing the Optimal Information Transport
on smooth, compact, and connected 2D manifolds $M$. We then provide the first
comparison and computations of solving the adaptive mesh methods on these
manifolds via Optimal Transport and Optimal Information Transport. In Section
1 we introduce the applied problem of adaptive moving mesh methods. These have
specific applications in the context of meteorology for the spherical
geometry, but have applications also in medical image registration and
computer graphics beyond the case of the sphere. In Section 2 we review
briefly the problem of Optimal Transport on the sphere (and beyond) and
highlight the regularity issues that arise which complicate the numerical
convergence guarantees. In Section 3 we outline how the Optimal Information
Transport PDE is derived. We also demonstrate that a reformulation of the
Poisson equation using the maximum operator will allow for the convergence
framework in [14] to apply to solving this PDE without boundary conditions. In
Section 4, we show the grid setup and how monotone operators are constructed
for the PDE in Optimal Transport and Optimal Information Transport on the
sphere. Then, we introduce the algorithm for the Optimal Information Transport
problem. In Section 5 we show computations for both problems on the sphere and
introduce how such computations can be generalized to smooth, compact, and
connected 2D manifolds $M$.

Acknowledgement: I would like to express my gratitude to Brittany Froese
Hamfeldt for fruitful discussions and for encouraging me to write this
manuscript on my own.

##  1\. Moving Mesh Methods

For this manuscript, our major application will be moving mesh methods, a type
of adaptive mesh method that has found special usage in computational PDE
techniques. Moving mesh methods require that the mapping be a diffeomorphism,
which means that the Jacobian of the mapping $T$ satisfies $0<J(T)<\infty$ in
the strong sense. This condition will prevent the mesh from tangling. The
resulting computed mesh is then often used in high-resolution PDE
computations, such as the Eady problem, which is a 2D vertical cross-section
of an incompressible Bousinnesq fluid, which was treated in [3]. Fixing the
number of grid points leads to simpler data structures (than some other
adaptive mesh techniques) in the PDE solving step, as was stressed in the
manuscripts [3, 31, 6].

In the Euclidean setup, the moving mesh problem setup requires the “physical”
target domain $\Omega_{p}\subset\mathbb{R}^{d}$ is where the PDE is posed with
coordinates $y$, while the “computational” input domain
$\Omega_{c}\subset\mathbb{R}^{d}$ with coordinates $x$ is usually chosen to be
a uniform rectangular grid (in Euclidean space). In non-Euclidean geometries,
one would like to use any off-the-shelf mesh generator for the computational
mesh. Then, one would like to find a diffeomorphic mapping from the off-the-
shelf “computational” mesh to the “physical” mesh $T:\Omega_{c}\to\Omega_{p}$.
Typically, in applications, the local density of grid points in the physical
domain $\Omega_{p}$ is determined by, for example, the time scales involved in
the solution of a fluid mechanics PDE. For example, in an evolving front or
shock, it is desirable for the density of points in $\Omega_{p}$ to be greater
than in the relatively unchanging parts of the solution, as in meteorological
applications where it is desirable to have high resolution in areas of high
precipitation [31]. This process could be done iteratively to get a very
accurate solution of the PDE, by solving the shock on more and more resolved
grids as one iterates. The way this is done is by feeding the information from
the PDE into a scalar monitor function $\mathcal{M}(y,t)>0$ and then solving

(1) $\mathcal{M}(y,t)J(T(x))=\theta(t)$

where $J$ is the Jacobian of the diffeomorphic mapping $T$ and
$\int_{\Omega_{p}}\mathcal{M}(y,t)dy=\theta(t)$. An example of a monitor
function from a function $f(y,t)$ is the scaled arc-length function:

(2) $\mathcal{M}(y,t)=\sqrt{1+\left|S\nabla_{y}f(y,t)\right|^{2}}$

where $S$ is a normalization factor [3]. The equation (1) is a change of
variables formula.

Applications of adaptive mesh techniques are also found in medical image
registration and sampling techniques, which require computing diffeomorphic
mapping between possibly very different density functions. Thus, we need to
solve the problem with a broad variety of applications in mind. The
diffeomorphic density matching problem can be solved via the Jacobian
equation. The existence of such a map is quite general, see [30].

######  Theorem 1 (Jacobian).

Let $M$ be an $n$-dimensional Riemannian manifold. Let $\mu_{0}$ and $\mu_{1}$
be two probability measures on $M$, and let $T:M\rightarrow M$ be a measurable
function such that $T_{\\#}\mu_{0}=\mu_{1}$. Let $\nu$ be a reference measure,
of the form $\nu(dx)=e^{-V(x)}\text{vol}(dx)$, where $V$ is continuous and vol
is the volume measure on $M$. Further assume that
$\mu_{0}(dx)=f_{0}(x)\nu(dx)$ and $\mu_{1}(dy)=f_{1}(y)\nu(dy)$, $T$ is
injective, and the distributional derivative of the mapping $DT$ is a locally
integrable function (this can be relaxed slightly). Then, $\mu_{0}$-almost
surely,

(3) $f_{0}(x)=f_{1}(T(x))J(T(x))$

where $J$ is the Jacobian determinant of $T$ at $x$, defined weakly by:

(4) $J(T(x)):=\lim_{\epsilon\rightarrow
0}\frac{\nu\left[T(B_{\epsilon}(x))\right]}{\nu[B_{\epsilon}(x)]}$

While the above theorem 1 stipulates the conditions for which we have a
distributional solution of the monitor equation (1), it does not state that
such a mapping $T$ is unique. We wish to find a particular mapping $T$ with a
computational scheme and have some a priori knowledge about it. In general, we
are not guaranteed a unique solution to (1). Hence, we need to impose further
conditions. A popular choice for regularity considerations is to require that
the map $T$ also minimizes the following integral:

(5) $I=\int_{\Omega_{c}}\left|T(x)-x\right|^{2}dx$

This choice leads to the mapping $T$ that arises in Optimal Transport, with
squared geodesic cost. It can be shown that such a mapping is irrotational
(meaning $\nabla\times T=0$ [4]) if the underlying manifold is
$\mathbb{R}^{d}$. Such an irrotational map is the gradient of a convex
function $u$, that is: $T(x)=\nabla_{x}u(x)$.

Injectivity of the map $T$ is guaranteed provided that the source and target
masses $\rho_{0}$ and $\rho_{1}$ are absolutely continuous with respect to
Lebesgue measure (or volume forms on manifolds), see [28]. Furthermore, the
source and target masses being absolutely continuous with respect to Lebesgue
measure ensures that the Optimal Transport problem is solved in the sense of
Aleksandrov, which means that $DT$ is locally integrable in the Euclidean case
[28]. In this case then, provided that $0<f_{0}<\infty$ and $0<f_{1}<\infty$,
we have $0<J(T)<\infty$, as can be seen by the equation (3). We still need
some regularity on the source and target probability measures in order to get
differentiability. The regularity theory of Caffarelli, see the summary in
[28], shows us that if $f_{0},f_{1}\in C^{0,\alpha}(\mathbb{R}^{d})$ for
$0<\alpha<1$, then we are guaranteed that the mapping is, in fact, a
diffeomorphism.

In the more general manifold setting, we have that the mapping solving the
Optimal Transport problem with squared geodesic cost is given by
$T(x)=\text{exp}_{x}\left(\nabla u\right)$ (as long as the manifold is
geodesically complete), see [22], where $u$ is known as a $c$-convex function,
a technical definition that is common in the Optimal Transport literature, see
[18] for a discussion of this condition and more about the general case of
manifolds. Thus, the question of whether the mapping has non-zero Jacobian
relates to the injectivity radius on the manifold $M$. For the $n$-sphere, the
injectivity radius is half the circumference of the sphere and fairly general
conditions on the source and target masses are given in [19], which show that
mass transport to a distance bounded strictly below the injectivity radius
gives us differentiability. In order to achieve the bounds $0<J(T)<\infty$,
again, we need to have $0<f_{0}<\infty$ and $0<f_{1}<\infty$. In the case of
the sphere, therefore, we can avoid tangling as well. More detail on the
regularity results of the Optimal Transport problem will be shown in 2. The
more general manifold case fails quite spectacularly in terms of
differentiability, see the result in [9] for a simple example where
differentiability fails to hold.

Optimal Information Transport provides another option for computing a mapping
for moving mesh methods. We will provide much more detail on the Optimal
Information Transport problem in Section 3, but we summarize the relevant
details here. The Optimal Information Transport formulation will be posed for
source and target probability distributions (volume forms) on a compact,
connected smooth manifold $M$ which have $C^{\infty}$ density functions that
are strictly bounded away from zero, and which leads to the existence of
diffeomorphic mappings $T$, which would prevent tangling.

##  2\. Optimal Transport on Compact Manifolds

In this section, we consider the Monge problem of Optimal Transport on $C^{3}$
compact, connected manifolds $M$ without boundary. That is, given source and
target probability densities $\mu_{0}$ and $\mu_{1}$ on the manifold that have
density functions $f_{0}$ and $f_{1}$, respectively, that is
$d\mu_{0}(x)=f_{0}(x)\text{dvol}(x)$ and $d\mu_{1}(x)=f_{1}(y)\text{dvol}(y)$
and a cost function $c(x,y):M\times M\rightarrow\mathbb{R}$, find a mapping
$T:M\rightarrow M$ such that $T_{\\#}\mu_{0}=\mu_{1}$, that is $\mu_{1}$ is
the pushforward measure of $\mu_{0}$ via the mapping $T$ such that:

(6) $T=\operatorname*{argmin}_{S}\int_{M}c(x,S)f_{0}(x)dx$

The Monge problem of Optimal Transport can be formulated as a partial
differential equation given some additional assumptions. We consider cost
functions $c(x,y)$ satisfying the following hypotheses 1, which are derived
from the hypotheses given in [21]:

######  Hypothesis 1 (Conditions on cost function).

We require the cost function satisfy the following conditions:

  1. (a)

The function $\nabla_{M,x}c(x,y)$ is Lipschitz for all $x,y\in M$ and
invertible in $y$

  2. (b)

For any $x,y\in M$, we have $\det D^{2}_{M,xy}c(x,y)\neq 0$

The optimal map corresponding to this cost function is determined from the
conditions

(7) $\begin{cases}\nabla_{M,x}c\left(x,T(x,p)\right)=-p,&x\in
M,p\in\mathcal{T}_{x}\\\ T(x,p)\in M\end{cases}\vspace*{-4pt}$

where $\mathcal{T}_{x}$ denotes the tangent plane at $x$. The solution to the
optimal transport problem is then given by

(8) $F\left(x,\nabla_{M}u(x),D_{M}^{2}u(x)\right)=0$

where

(9) $F(x,p,M)\equiv-\det\left(M+A(x,p)\right)+H(x,p)$

subject to the $c$-convexity condition, which requires

(10) $D_{M}^{2}u(x)+A(x,\nabla_{M}u(x))\geq 0.$

Here

$\displaystyle A(x,p)$ $\displaystyle=D_{M,xx}^{2}c\left(x,T(x,p)\right)$
$\displaystyle H(x,p)$
$\displaystyle=\left|\det{D_{M,xy}^{2}c\left(x,T(x,p)\right)}\right|f_{0}(x)/f_{1}\left(T(x,p)\right),$

and the PDE now describes a nonlinear relationship between the surface
gradient and Hessian on the manifold.

###  2.1. Regularity

We consider the optimal transport problem (8) under the following hypothesis.

######  Hypothesis 2 (Conditions on data (smooth)).

We require problem data to satisfy the following conditions:

  1. (a)

There exists some $m>0$ such that $f_{1}(x)\geq m$ for all
$x\in\mathbb{S}^{n}$.

  2. (b)

The mass balance condition holds,
$\int_{\mathbb{S}^{2}}f_{0}(x)\,dx=\int_{\mathbb{S}^{2}}f_{1}(y)\,dy$.

  3. (d)

The data satisfies the regularity requirements $f_{0},f_{1}\in
C^{1,1}(\mathbb{S}^{n})$.

  4. (e)

The cost function $c(x,y)$ is either the squared geodesic cost
$c(x,y)=\frac{1}{2}d(x,y)^{2}$, the logarithmic cost $c(x,y)=-2\log(1-x\cdot
y)$, or $f(d(x,y))$ for some $f:[0,\pi)\rightarrow\mathbb{R}$ smooth and
strictly increasing with $f^{\prime}(0)=0$

  5. (f)

The geometry is the $n$-sphere $\mathbb{S}^{n}$

We thus get the following regularity result, from Loeper [19]

######  Theorem 3 (Regularity).

The optimal transport problem (8) with data satisfying Hypothesis 2 has a
$c$-convex classical solution $u\in C^{3}(\mathbb{S}^{2})$.

The solution to (8) is unique only up to additive constants. In this
manuscript, sometimes we will fix a point $x_{0}\in M$ and add the additional
constraint:

(11) $u(x_{0})=0.$

At other times, it may be more convenient to choose the mean-zero solution,
that is to impose the constraint

(12) $\int_{M}u=0.$

###  2.2. Beyond the Sphere

In some cases, manifolds $M$ with non-positive cost-sectional curvature at any
point $x\in M$ can have positive measures $\mu_{0},\mu_{1}\in C^{\infty}(M)$,
but $T$ is not even guaranteed to be continuous [18]. In order to explain what
the cost-sectional curvature is, we introduce the Ma, Trudinger, Wang tensor
[21, 18]:

(13)
$\mathfrak{G}_{c}(x_{0},y_{0})(\xi,\nu)=D^{4}_{p_{\nu},p_{\nu},x_{\xi},x_{\xi}}\left[(x,p)\mapsto-c\left(x,T_{x_{0}}(p)\right)\right]|_{x_{0},p_{0}=-\nabla_{x}x(x_{0},y_{0})}$

The cost-sectional curvature is negative at a point $(x,y)$ if there exist
$\xi,\nu$ such that
$\mathfrak{G}_{c}(x,y)(\xi,\nu)<C_{0}\left|\xi\right|^{2}\left|\nu\right|^{2}$.
The idea is more perhaps transparent if one looks at the squared geodesic cost
case. In this simple case, the negativity of the cost-sectional curvature is
equivalent to the negativity of the sectional curvature of $M$ due to the
equality:

(14)
$\frac{\mathfrak{G}_{c}(x,x)(\nu,\xi)}{\left|\xi\right|^{2}\left|\nu\right|^{2}-(\xi\cdot\nu)^{2}}=\frac{2}{3}\cdot\text{sectional
curvature of $M$ at $x$ in the plane $(\xi,\nu)$}$

But even if the underlying manifold has strictly positive curvature
everywhere, the Ma, Trudinger, Wang tensor $\mathfrak{G}_{c}(x,y)$ may be non-
negative on the off-diagonal $(x,y)\ \text{s.t}\ x\neq y$. This was shown true
even for some ellipsoids of revolution in [9]. In some cases (such as $C^{4}$
perturbations from the $n$-sphere [9]) we still can have classical solutions
given the correct hypotheses on the cost function and source and target mass
density functions.

##  3\. Optimal Information Transport on Compact Manifolds

The Monge problem of Optimal Transport hints at a deeper relation between the
space of probability distributions and the mappings between probability
measures. We have seen in section 2 that the unique mapping $T$ of the Monge
problem is used to define the total cost of transporting the probability
measure $\mu_{0}$ to $\mu_{1}$. If the cost function is chosen to be the
squared geodesic cost, this total cost defines a distance between probability
measures, see [28, 30], usually referred to as the Wasserstein distance.
Furthermore, a Wasserstein interpolation between the source mass $\mu_{0}$ and
the target mass $\mu_{1}$, denoted by $\mu(t)$ can be uniquely defined from a
convex combination of the identity map and the Optimal Transport map from
$\mu_{0}$ to $\mu_{1}$, that is: $T(t)=(1-t)\text{Id}+tT$, by defining
$\mu(t)=T(t)_{\\#}\mu_{0}$, which defines a geodesic in Wasserstein space.

Another alternative distance between probability measures arising primarily in
statistical applications is the Fisher-Rao distance. This distance is defined
as the Riemannian distance arising from the Fisher-Rao metric on the space of
probability measures. The Fisher-Rao metric, in turn, is the Hessian of the
Kullback-Liebler divergence between the probability measures $\mu$ and $\nu$,
that is the quantity

(15) $\text{KL}(\mu,\nu):=\int_{M}\log\left(\frac{d\mu}{d\nu}\right)d\mu$

which measures the relative entropy between one probability measure $\mu$ and
another $\nu$. Given the space of probability measures then equipped with the
Fisher-Rao metric, one may ask the question if there is a mapping $T$ which
characterizes the Fisher-Rao distance in the same way that the Optimal
Transport mapping $T$ characterizes the Wasserstein distance and
interpolation. The short answer is yes, and one can derive it by analyzing the
geometric connection between the space of diffeomorphisms $\text{Diff}(M)$ and
the space of positive smooth probability measures (volume forms)
$\text{Dens}(M)$ on the underlying manifold $M$. That is, one can define a
projection map $\pi:\text{Diff}(M)\rightarrow\text{Dens}(M)$ providing a
principle bundle structure, introduced by Moser [25]. This structure is then
analyzed to show that, sure enough, a unique Optimal Information Transport map
$T$ can be associated with the Fisher-Rao distance between $\mu_{0}$ and
$\mu_{1}$ and define a unique interpolation as well, and a simple system of
PDE that can solved for $T$.

For more information about Optimal Information, see the background and more
detail presented in [1, 23]. In following section, we present a brief summary
of the derivation of the Optimal Information Transport PDE. For all the
discussion in this section, we will assume that the underlying manifold $M$ is
compact and connected and is equipped with the standard volume form vol. We
introduce the Fréchet Manifold of smooth volume forms over $M$ with total
volume $\text{vol}(M)$:

(16)
$\text{Dens}(M)=\left\\{\mu\in\Omega^{n}(M):\int_{M}\mu=\text{vol}(M),\mu>0\right\\}$

where the notation $\Omega^{n}(M)$ denotes the $n$-forms over the manifold
$M$, using topology induced by the Sobolev seminorms. Although the Fisher-Rao
metric:

(17)
$G^{F}_{\mu}(\alpha,\beta)=\frac{1}{4}\int_{M}\frac{\alpha}{\mu}\frac{\beta}{\mu}\mu$

is weak in the Fréchet topology, the geodesics are explicit!

We introduce another manifold with rich structure. The set of diffeomorphisms
on $M$ is denoted $\text{Diff}(M)$ with topology induced by the Sobolev
seminorms. We have assumed $M$ compact, which makes the set $\text{Diff}(M)$
is a Fréchet Lie group under the composition of maps. The Lie algebra of
$\text{Diff}(M)$ is given by the space $\mathfrak{X}(M)$ of smooth vector
fields. We equip this manifold with the information metric on
$\text{Diff}(M)$, defined by:

(18) $G^{I}_{\varphi}(U,V)=\int_{M}g(\Delta
u,v)\text{vol}+\lambda\sum_{i=1}^{k}\int_{M}g(u,\xi_{i})\text{vol}\int_{M}g(v,\xi_{i})\text{vol}$

This is also weak in the Fréchet topology, but its geodesics are well-posed.

Now, the diffeomorphic maps define a group action on the space of densities,
via the pullback $\varphi^{\\#}\mu=\nu$. The volume preserving maps (with
respect to $\varphi$): define the isotropy group:

(19)
$\text{Diff}_{\mu}(M):=\left\\{\varphi\in\text{Diff}(M);\varphi^{\\#}\mu=\mu\right\\}$

This group action is transitive, that is, given $\mu$ and $\nu$, there exists
a diffeomorphism $\varphi$ such that $\varphi^{\\#}\mu=\nu$. Therefore, the
manifold of densities $\text{Dens}(M)$ is isomorphic to the space of right
cosets $\text{Diff}(M)/\text{Diff}_{\mu}(M)$, and hence the density
$\mu\in\text{Diff}(M)$ parametrizes a partition of the space of densities.
Furthermore, by defining the projection map:

(20)
$\pi_{\mu}:\text{Diff}(M)\ni\varphi\mapsto\varphi^{\\#}\mu\in\text{Dens}(M)$

we can define the following principal bundle structure:

${\text{Diff}_{\mu}(M)}$${\text{Diff}(M)}$${\text{Dens}(M)}$$\scriptstyle{\pi_{\mu}}$

The goal is now to give a Riemannian structure to the principal fiber bundle
using (17) and (18). The upshot is that although $\pi$ is a submersion it is
actually also a Riemannian submersion with respect to $G^{F}$ and $G^{I}$,
i.e. the metric $G^{I}_{\varphi}$ descends to $G^{F}$. That is:

(21)
$G^{I}_{\varphi}(U,V)=G^{F}_{\pi_{\mu}(\varphi)}\left(T_{\varphi}\pi_{\mu}\cdot
U,T_{\varphi}\pi_{\mu}\cdot V\right)$

Descending metrics “descend” horizontal geodesics, that is if $\varphi(t)$ is
a horizontal geodesic curve in $\text{Diff}(M)$, then
$\mu(t):=\pi(\varphi(t))$ is a geodesic curve in $\text{Dens}(M)$.

The problem of diffeomorphic density matching is that given a path of
densities $\mu(t)$, we desire to find the horizontal path of diffeomorphisms
$\varphi(t)$ which project onto $\mu(t)$, that are also of minimal length with
respect to $G^{I}$. That is, solve the exact density matching problem, that is
find a $\varphi(t)$ such that:

(22) $\begin{cases}\varphi(0)=\text{id},\\\ \varphi^{\\#}\mu_{0}=\mu(t)\\\
\text{minimizing}\
\int_{0}^{1}G^{I}_{\varphi(t)}\left(\dot{\varphi}(t),\dot{\varphi}(t)\right)dt\end{cases}$

From now on, we assume that $\text{vol}=\mu_{0}$, i.e. $\mu_{0}$ is the
canonical volume form. We take the equation
$\varphi^{\\#}(t)\mu_{0}=\varphi^{\\#}(t)\text{vol}=\mu(t)$ and differentiate
with respect to $t$:

(23)
$\dot{\mu}(t)=\partial_{t}\left(\varphi(t)^{\\#}\text{vol}\right)=\varphi^{\\#}\text{div}_{\text{vol}}v(t)$

where $v(t)=\dot{\varphi}\circ\varphi^{-1}$. This can be rewritten as:

(24) $\dot{\mu}(t)=\text{div}\left(v(t)\right)\circ\varphi(t)\mu(t)$

From here we perform the Hodge-Helmholz decomposition for the vector field
$v$, by writing $v=\text{grad}f+w$. It turns out that the Hodge-Helmholtz
composition is orthogonal with respect to the information metric $G^{I}$ so
the length of the path $\varphi(t)$ is minimal for $w=0$. Therefore, we must
solve the following Poisson equation for the curl-free term $f$:

(25) $\begin{cases}\Delta
f(t)=\frac{\dot{\mu}(t)}{\mu(t)}\circ\varphi(t)^{-1}\\\
\dot{\varphi}(t)=\text{grad}\left(f(t)\right)\circ\varphi(t),\quad\varphi(0)=\text{id}\end{cases}$

As noted before, the geodesics for the Fisher-Rao metric have explicit forms!
Thus, we can insert explicit forms for $\mu(t)$ and $\dot{\mu(t)}$. Define
$W:\text{Dens}(M)\rightarrow C^{\infty}(M)$ by
$\mu\mapsto\sqrt{\frac{\mu}{\text{vol}}}$, then the geodesics are given by:

(26) $[0,1]\ni
t\mapsto\left(\frac{\sin\left((1-t)\theta\right)}{\sin\theta}f_{0}+\frac{\sin(t\theta)}{\sin\theta}f_{1}\right)^{2}\text{vol}$

where

(27) $\theta=\arccos\left(\frac{\left\langle
f_{0},f_{1}\right\rangle_{L^{2}}}{\text{vol}(M)}\right)$

and $f_{i}=W(\mu_{i})$. The Optimal Information Transport problem is then
given by

(28) $\begin{cases}\varphi(0)=\text{id},\\\
\varphi_{\\#}\mu_{0}=\mu(t)\end{cases}$

where $\varphi_{\\#}\mu_{0}$ denotes the pushforward of $\mu_{0}$. The
solution of this problem (28) is then given by $\varphi^{-1}(1)$, where
$\varphi(1)$ solves Equation (25). We also have a solution for the inexact
problem, where the object is to find $\varphi$ that minimizes:

(29) $\sigma
d^{2}(\text{id},\varphi)+\bar{d}^{2}(\varphi^{\\#}\mu_{0},\mu_{1})^{2},\quad\sigma>0$

where

(30)
$d^{2}(\varphi,\psi):=\inf_{\gamma(0)=\varphi,\gamma(1)=\psi}\int_{0}^{1}G_{\gamma(t)}\left(\dot{\gamma}(t),\dot{\gamma}(t)\right)dt$

and $\bar{d}$ denotes the distance function on $\text{Dens}(M)$ with
Riemannian metric $\bar{G}$. Thus, the first term in (30) is a regularity term
and the second term penalizes $\varphi$ if the pushforward of $\mu_{0}$ is too
far from $\mu_{1}$. Let $\varphi(1)^{-1}$ solve the exact problem (22) then
let $s=\frac{1}{1+\sigma}$. The solution to the inexact problem is then given
by $\varphi(s)^{-1}$.

##  4\. Numerical Methods for Optimal Transport and Optimal Information
Transport

###  4.1. Convergent Numerical Methods for Elliptic PDE on Compact Manifolds

We assume that we are on a compact, connected two-dimensional $C^{3}$ manifold
$M$. The convergence framework and schemes we describe here are inspired by a
numerical scheme recently designed by the authors for the $2$-sphere, which is
equipped with a proof of convergence to the physically meaningful solution of
the Optimal Transport problem, see [14, 15] for details and analysis. The
adaptation of the convergence framework from the sphere to $M$ requires little
modification provided that the solution to the elliptic PDE in question is
merely a priori know to be at least Lipschitz.

The discretizations presented here must be consistent and monotone. We also
get compactness from Lipschitz stability. We will have in mind that we are
using finite-difference schemes, but actually the convergence framework is
more general than this. For the sphere, in the manuscript [15] finite-
difference schemes were, in fact, constructed which satisfied the hypotheses
of the convergence framework. We have in mind that, similarly, finite-
difference schemes can be constructed which satisfy the hypotheses for more
general manifolds $M$, but defer the explicit construction of such schemes to
later work. The properties of consistency, monotonicity, and Lipschitz
stability are sufficient to guarantee convergence in nonsmooth cases and can
be used to build explicit convergence rates, an avenue of research that is
being pursued in a more general case (for the Optimal Transport PDE) in
another manuscript by the author.

The convergence framework hinges first on the reformulation of the PDE to
include some control on its Lipschitz constant. Given the underlying elliptic
PDE that we wish to solve (such as a Monge-Ampére-type PDE for the Optimal
Transport problem or a Poisson PDE for the Optimal Information Transport
problem) denoted by

(31) $F(x,Du(x),D^{2}u(x))=0$

we instead discretize the following modified PDE:

(32) $\max\\{F(x,\nabla u(x),D^{2}u(x),\left\|\nabla u\right\|-R\\}=0$

where $R$ is chosen to be strictly greater than the a priori Lipschitz bound.
After discretization, the Lipschitz bound imparts sufficient stability on the
solution to prove overall uniform convergence of the discrete solution $u^{h}$
to the actual solution $u$, provided that we have constructed a consistent,
monotone scheme. The equivalence of the solution of (32) to (31) for the
Optimal Transport case was presented in the manuscript [14]. Here we present
the argument for the Poisson equation.

First, we introduce the notion of a continuous viscosity solution, since the
compactness arguments are very clean. The arguments for uniform convergence
for elliptic PDE also are more naturally extended to weaker cases using the
notion of viscosity solutions. Suppose we have the elliptic PDE operator

######  Definition 4 (Continuous Viscosity Solutions).

An upper (lower) semicontinuous function $u:M\rightarrow\mathbb{R}$ is a
viscosity sub-(super)- solution of the PDE (31) if for every $x_{0}\in M$ and
$\phi\in C^{\infty}(M)$ such that $u-\phi$ has a local maximum (minimum) at
$x_{0}$ we have:

(33) $F\left(x_{0},u(x_{0}),\nabla u(x_{0}),D^{2}u(x_{0})\right)\leq(\geq)0$

A continuous function $u:M\rightarrow\mathbb{R}$ is a viscosity solution of
(31) if it is both a sub- and super-solution.

The notion of test functions is slightly different in the case of the Optimal
Transport PDE, as defined in [14]. Note also that, critically, the condition
(33) can hold vacuously. We modify the Poisson equation by defining the
operator

(34) $F(x,u(x),\nabla u(x),D^{2}u(x)):=\max\left\\{-\Delta
u+f(x),\left\|u(x)\right\|-R\right\\}=0$

The function $f(x)$ satisfies the compatibility condition: $\int_{M}f(x)dx=0$.
We now state the equivalence theorem that shows that adding a term to the
modified Poisson equation (34) does not change the solution. The proof of this
theorem is given in Appendix Appendix A: Modified Poisson.

######  Theorem 2.

For any sequence $\epsilon_{n}\rightarrow 0$, $\epsilon_{n}>0$ and for $R>0$
large enough, the Lipschitz viscosity solution $w_{\epsilon_{n}}$

(35) $\max\left\\{-\Delta
w_{\epsilon_{n}}(x)+f(x),\left\|w_{\epsilon_{n}}(x)\right\|-R\right\\}+\epsilon_{n}w_{\epsilon_{n}}(x)=0,$

where $f(x)$ satisfies $\int_{M}f(x)dx=0$, converges to the mean-zero
Lipschitz viscosity solution $u$ of the Poisson equation:

(36) $-\Delta u(x)+f(x)=0$

This now allows us to instead discretize the PDE (34) and likewise for the
Optimal Transport problem.

Now, we describe the grid and discretization on the manifold. We start with a
finite set of $N$ grid points $\mathcal{G}\subset M$ that lie on the manifold.
We let $d_{M}(x,y)$ denote the usual geodesic distance between points $x,y$ on
the manifold $M$. To the grid $\mathcal{G}$, we associate a number $h$ which
serves as a parameter indicating the overall spacing of grid points. More
precisely,

(37) $h=\sup\limits_{x\in
M}\min\limits_{y\in\mathcal{G}^{h}}d_{M}(x,y)=\mathcal{O}\left(N^{-1/2}\right).$

In particular, this guarantees that any ball of radius $h$ on the sphere will
contain at least one discretization point. The grid must sufficiently resolve
the manifold and satisfy some technical regularity requirements which are not
difficult in practice to achieve, as detailed in [14]. These regularity
requirements ensure that the grid is more or less isotropic with respect to
direction. Now we seek an approximation of the form

(38) $F^{h}(x,u;f_{1},f_{2})=0,\quad x\in\mathcal{G}$

At this point, we will begin to more explicitly state the construction of the
finite-difference scheme in the case of the sphere. We defer its full
generalization to compact surfaces $M$ to future work. We now consider a fixed
grid point $x_{i}\in\mathcal{G}$ and a grid function
$u:\mathcal{G}\to\mathbb{R}$ and explain how we obtain the value of
$F^{h}(x_{i},u;f_{1},f_{2})$. To perform the computation at $x_{i}$, we begin
by projecting grid points close to $x_{i}$ onto the tangent plane at $x_{i}$.
That is, we consider the set of relevant discretization points

(39) $\mathcal{Z}(x_{i})=\left\\{z=\text{Proj}(x;x_{i})\mid
x\in\mathcal{G},d_{M}(x,x_{i})\leq\sqrt{h}\right\\}.$

In order to minimize the consistency error of the overall scheme, the choice
of the radius of the neighborhood is $\sqrt{h}$ is optimal, see [10], which
will allow the second directional derivative to eventually resolve every
direction in the computational neighborhood.

The projection onto local tangent planes is accomplished using geodesic normal
coordinates, which are chosen to preserve the distance from $x_{i}$ (i.e.
$d_{M}(x,x_{i})=\|x_{i}-\text{Proj}(x;x_{i})\|$). This prevents any
distortions that the non-zero curvature of a manifold necessarily introduces
to second-order derivatives. This allows for monotone approximation schemes to
be more easily constructed at least in the case of the sphere. Beyond the case
of the sphere, it is necessary to merely have a consistent approximation for
local geodesic coordinates, since an explicit formula does not exist in
general.

The Optimal Transport PDE requires computing gradients and eigenvalues of
Hessian matrices. In order to compute the gradient, for example, we need to
search various directions $\nu$ in the neighborhood. For the Optimal
Information Transport problem we also need to compute a gradient. Therefore,
we will consider the following finite set of possible directions,

(40)
$V=\left\\{\left\\{(\cos(jd\theta),\sin(jd\theta)),(-\sin(jd\theta),\cos(jd\theta))\right\\}\mid
j=1,\ldots,\frac{\pi}{2d\theta}\right\\},$

where the angular resolution
$d\theta=\dfrac{\pi}{2\lfloor\pi/(2\sqrt{h})\rfloor}$.

In order to build approximations for derivatives in each direction $\nu\in V$,
we will select four grid points $x_{j}\in\mathcal{Z}(x_{i})$, $j=1,\ldots,4$,
which will be used to construct the directional derivatives in this direction.
To accomplish this, we let $\nu^{\perp}$ be a unit vector orthogonal to $\nu$
and represent points in $x\in\mathcal{Z}(x_{i})$ using (rotated) polar
coordinates $(r,\theta)$ centred at $x_{i}$ via

$x=x_{i}+r(\nu\cos\theta+\nu^{\perp}\sin\theta),x\in\mathcal{Z}(x_{i}).$

Then we select four points, each in a different quadrant
($Q_{1},\ldots,Q_{4}$), that are well-aligned with the direction of $\nu$ and
approximately symmetric about the axis $\nu$ via

(41)
$x_{j}\in\operatorname*{argmin}\limits_{x\in\mathcal{Z}(x_{i})}\left\\{\left|\sin\theta\right|\mid\left|\sin\theta\right|\geq
d\theta,r\geq\sqrt{h}-2h,x\in Q_{j}\right\\}$

where $\cos\theta\geq 0$ for points in $Q_{1}$ or $Q_{4}$ and $\sin\theta\geq
0$ for points in $Q_{1}$ or $Q_{2}$, see Figure 1.

Figure 1. The computational points are chosen so as to ensure that they are
within $2d\theta$ of the axis of computation, but also stay at least $d\theta$
away from the axis to ensure monotonicity of the scheme

From here, we construct approximations of second directional derivatives (and
first directional derivatives for the usual coordinate directions $(1,0)$ and
$(0,1)$) of the form

(42)
$\begin{split}\mathcal{D}_{\nu\nu}u(x_{i})&=\sum\limits_{j=1}^{4}a_{j}(u(x_{j})-u(x_{i}))\approx\frac{\partial^{2}u(x_{i})}{\partial\nu^{2}}\\\
\mathcal{D}_{\nu}u(x_{i})&=\sum\limits_{j=1}^{4}b_{j}(u(x_{j})-u(x_{i}))\approx\frac{\partial
u(x_{i})}{\partial\nu}.\end{split}$

After performing a Taylor expansion of the sums in (42) we get explicit values
for $a_{j}$ and $b_{j}$, which can be found in [15, 16]. The Laplacian can be
built from approximations of second-order directional derivatives. Because the
Laplacian is the trace of the Hessian, we take a geodesic normal coordinate
system $(x^{j})_{j}$ and define:

(43)
$\Delta^{h}u(x_{i})=\mathcal{D}_{\nu_{x_{1}}\nu_{x_{1}}}u(x_{i})+\mathcal{D}_{\nu_{x_{2}}\nu_{x_{2}}}u(x_{i})$

In the case where the underlying PDE contains functions of the gradient, the
goal is to add or subtract a Laplacian regularizing term in order to build
monotonicity into the discretization, more details of which can be found in
[15]. The full discretization for the discrete operator $F^{h}$ for the
Optimal Transport problem on the sphere then is:

(44) $F^{h}\left(x_{i},u^{h}(x_{i})\right)=\\\ \min_{\\{\nu_{1},\nu_{2}\\}\in
V}\prod_{j=1}^{2}\max\left\\{\mathcal{D}^{h}_{\nu_{j}\nu_{j}}u^{h}(x_{i})+g^{-}_{1,\nu_{j}}\left(\nabla^{h}u^{h}(x_{i})\right),0\right\\}-f_{1}(x_{0})g^{+}_{2}\left(\nabla^{h}u^{h}(x_{i})\right)$

where $g_{1}(p;\nu)=-\mathcal{D}^{h}_{\nu\nu}c(x_{i},y)|_{y=T(x_{i},p)}$ and
$g_{2}(p;\nu)=\frac{\sin\left\|p\right\|}{\left\|p\right\|}$ and
$g^{\pm}_{i,\nu}\left(\nabla^{h}u^{h}(x_{i})\right)=g_{i}\left(\nabla^{h}u^{h}(x_{i});\nu\right)\mp\epsilon_{g}\Delta^{h}u^{h}(x_{i})$

And the full discretization for the discrete operator $F^{h}$ solving a
Poisson equation on the sphere is:

(45)
$F^{h}\left(x_{i},u^{h}(x_{i})\right)=\Delta^{h}u^{h}(x_{i})+\epsilon^{h}u^{h}(x_{i})=f(x_{i})$

where $\epsilon^{h}\rightarrow 0$ is chosen to be the consistency error of the
operator $\Delta^{h}$. The discretization of the Lipschitz constraint is:

(46)
$E^{h}\left(x_{i},u^{h}(x_{i})\right):=\max\left|\nabla^{h}u^{h}(x_{i})\right|-R$

and we define:

(47)
$G^{h}\left(x_{i},u^{h}_{n}(x_{i})\right):=\max\left\\{F^{h}\left(x_{i},u^{h}_{n}(x_{i})\right),E^{h}\left(x_{i},u^{h}_{n}(x_{i})\right)\right\\}$

Then, the procedure for solving the problem is:

Algorithm 1 Computing the solution to elliptic PDE $F[u]=0$

1:Initialize $u^{h}_{0}$;

2:Fix $\epsilon>0$;

3:while $\left|G^{h}\left(x_{i},u^{h}_{n}(x_{i})\right)\right|>\epsilon$ do

4: Compute $u^{h}_{n+1}(x_{i})=u^{h}_{n}(x_{i})+\Delta
tG^{h}\left(x_{i},u^{h}_{n}(x_{i})\right)$

5:end while

######  Remark 5.

The computational method of [15] in principal involves solving the problem
with this approximation, verifying that the solution satisfies required
Lipschitz bounds, then if necessary solving a second discrete problem to
enforce the Lipschitz condition. However, we do not see the verification step
fail in practice, and hence never actually need to solve using the full
operator $G^{h}$, but instead simply use the operator $F^{h}$.

As can be seen from the algorithm, parabolic solvers are employed to solve the
discretized system, especially when $F^{h}$ is a nonlinear operator as was
shown in [15]. In the case of the Poisson equation when the Lipschitz
constraint is inactive, we solve the linear system via standard linear algebra
techniques.

###  4.2. Uniform Convergence Framework for the Poisson Equation for Compact
Surfaces without Boundary

The key properties in the discretization of the elliptic PDE without boundary
for both the Optimal Transport and Optimal Information Transport PDE are
consistency, monotonicity, and Lipschitz stability. Now, we show that with
these properties we can derive a uniform convergence result for the Poisson
equation on a compact 2D manifold without boundary. We use the auxiliary
numerical scheme to prove boundedness:

(48)
$\tilde{G}^{h}(v^{h})(x_{i}):=\max\left\\{-\Delta^{h}v^{h}(x_{i})-h^{\beta}v^{h}(x_{i})+f(x_{i}),E^{h}\left(x_{i},v^{h}(x_{i})\right)\right\\}+\\\
h^{\alpha}v^{h}(x_{i})=0$

Let the formal consistency error of the scheme be
$\mathcal{O}\left(h^{\alpha}\right)$. Fix a point $x_{0}\in M$. Choose the
solution $u(x)$ such that $u(x_{0})=0$ and define $u^{h}$ such that
$u^{h}(x_{0})=0$ by effecting the choice $u^{h}(x)=v^{h}(x)-v^{h}(x_{0})$.

######  Lemma 6 (Boundedness of Solution to Auxiliary Numerical Scheme).

The solution $v^{h}$ is bounded if $\beta\leq\alpha$.

###### Proof.

The discretization $\tilde{G}^{h}$ is monotone and proper. Therefore, it has a
discrete comparison principle, from the results in [26]. Then, we plug in
$v^{h}$ and $u\pm C$ and apply the discrete comparison principle. The uniform
bound follows from an identical argument for both the smooth and nonsmooth
cases as shown in [14] for the Optimal Transport PDE. ∎

From the details of the proof in [14] we choose $\beta\leq\alpha$. It turns
out that it is best to choose $\beta=\alpha$. We assume this choice from now
on. We also assume that we can talk about $u^{h}(x)$, for any $x\in M$,
assuming that we have used the interpolation of $u^{h}$ on the sphere detailed
in [14], with modifications, mutatis mutandis, to apply the interpolation to
the manifold $M$.

The uniform convergence critically depends on the solution of the modified PDE
(34) being the same as the solution of the Poisson equation. This was shown in
4, while a similar result for the Optimal Transport PDE is in the paper [14].
Thus, we have:

######  Theorem 3 (Uniform Convergence in Nonsmooth Cases).

Assuming that we have a consistent and monotone discretization of $\Delta^{h}$
and the Lipschitz control on the discretization, then the solution $u^{h}$
converges uniformly to the Lipschitz viscosity solution $u$ of the Poisson
equation on the manifold $M$ where $u(x_{0})=0$

###### Proof.

Techniques are identical to those used in [14]. Uniform convergence can be
proved for nonsmooth cases ($C^{0,1}(M)$) by assuming that the modified
Poisson equation (34) has a unique solution. ∎

######  Remark 7.

The compactness arguments given in this section naturally do not lead to
explicit convergence rates. However, forthcoming work will establish the
guaranteed convergence rate
$\left|u^{h}-u\right|=\mathcal{O}\left(h^{\alpha/3}\right)$ for a priori
$C^{3}$ solutions of the Poisson equation by explicitly constructing bounding
functions.

###  4.3. Algorithm for Optimal Information Transport

In order to solve the Optimal Information Transport problem on $M$, we have to
perform some pre-computations involving the geometry and the mass densities.
Then, we iterate an Euler scheme where each step involves one iteration of a
Poisson solver. Here are the steps to compute the approximation to the
diffeomorphic mapping $\varphi(1)$. Assume $\mu_{1}\in L^{\infty}(M)$. The
forward map at a time step $n$ will be denoted by $T_{n}$ and the
corresponding inverse map at a time step $n$ will be denoted by $S_{n}$. We
will outline the algorithm for the sphere. More general manifolds require some
fine-tuning, introduced in 5. The algorithm requires two important functions:
Proj a projection map consistent with the exponential map and Interp, a
consistent interpolation map. For the case of the sphere, with an explicit
formula for the exponential map, the function Proj is simple to execute.
However, the interpolation function is perhaps more challenging. A robust
interpolation is that described in [14], where the values are consistently
interpolated over the Delaunay triangles.

Algorithm 2 Computing the diffeomorphic mapping $\varphi=S(1)$

1:Initialize $T_{0}=\text{id}$;

2:Initialize $S_{0}=\text{id}$;

3:Fix $\Delta t\ll 1$;

4:Precompute $\theta$ via quadrature;

5:while $n\Delta t<1$ do

6: Compute the density function $\nu_{n}:=\dot{\mu}_{n}/\mu_{n}$ using the
explicit formulas;

7: Interpolate: $\text{Interp}\left\\{\nu_{n}\right\\}$ onto the grid
$\left\\{S(x_{j})\right\\}_{j}$;

8: Solve $\Delta^{h}f_{n}(x_{i})=\nu_{n}\left(S_{n}(x_{i})\right)$ with any
monotone, consistent discretization of (34);

9: Compute $\nabla^{h}f_{n}(x_{i})$ for all $x_{i}$;

10: Interpolate: $\nabla^{h}f_{n}(x_{i})$ onto the grid
$\left\\{T(x_{j})\right\\}_{j}$;

11: Compute $T_{n+1}(x_{i})=\text{Proj}\left\\{T_{n}(x_{i})+\Delta
t\nabla^{h}f_{n}\left(T_{n}(x_{i})\right)\right\\}$;

12: Compute $S_{n+1}(x_{i})=S_{n}\left(\text{Proj}\left\\{x-\Delta
t\nabla^{h}f_{n}(x_{i})\right\\}\right)$

13:end while

##  5\. Implementation

###  5.1. Improved Mesh Method Computations on the Sphere

Our objective is to perform adaptive mesh computations on the sphere. The
original mesh will be assumed to be given by some off-the-shelf mesh generator
which is simple to implement. Here we will be using a mesh defined on a cube
which is then projected onto the sphere. The original mesh, therefore, does
not have constant density. However, in the moving mesh problem, we will be
performing computations that produce diffeomorphic maps from a constant source
density to a variable target density. This allows us to redistribute the mesh
as desired.

First, we demonstrate that the computation of mesh redistribution using
Optimal Transport and Optimal Information Transport can be both successful in
produce meshes which do not exhibit tangling. We select density functions with
the goal of producing a transport map $T$ that will concentrate mesh points
around the equator, see (49).

(49) $\begin{cases}f_{1}(\theta,\phi)=\frac{1}{4\pi}\\\
f_{2}(\theta,\phi)=\left(1-\text{exp}\left(-1/30\left(\text{arccos}(z)-\pi/2\right)^{2}\right)\right)/3.53552\end{cases}$

The original mesh is generated from the cube mesh that is then projected onto
the sphere, see 2

Figure 2. A $N=5048$-point cube mesh from the top 2 and diagonally above
looking down 2

We then require that the density concentrates about the equator according to
(49), see the visual representation in the figure 3

Figure 3. Density in the original mesh (left) and the modified density (right)

The resulting mesh restructuring computed using Optimal Transport is pictured
in Figure 4.

Figure 4.  4 3D view of the mesh obtained via Optimal Transport map, 4 top
view, and a detailed view of the equator 4 showing that the grid lines do not
tangle

Now we perform the same computation using Optimal Information Transport, see
5:

Figure 5.  5 3D view of the mesh obtained via Optimal Information Transport, 5
top view, and a detailed view of the equator 5 showing that the grid lines do
not tangle

Another example serves to demonstrate the versatility and desirability of
using Optimal Information Transport for moving mesh problems. We take a
constant source density $f_{1}=1/(4\pi)$ and map it to a complicated,
discontinuous map of the world, where a higher density of grid points is
required over the oceans and lower density over the continents, see Figure 6.

Figure 6. Constant density in the source density (right) and the target globe
density (left)

As before, we use the cube mesh 2 and transform it via Optimal Information
Transport. The result is shown in Figure 7.

Figure 7. Entire spherical mesh transformed by Optimal Information Transport 7
in 3D view. The view is of Africa, Europe and Asia, but the other side of the
sphere mesh is also visible. Shown on the right is detail on North America
with other side of sphere hidden from view in order to show the mesh more
clearly 7. The outlines of North and South America can be discerned and
tangling is minimized.

Some difficulties are encountered when attempting to perform the same
numerical computation via Optimal Transport. First, the necessity of using an
interpolation function for the target mass density $g$ onto the mesh defined
by $T(x)$ leads to instabilities in the scheme due to the very nonlinear
nature of the equation we are trying to solve. Compensating for this, by using
an inverted image for the source mass $f$ and a constant for the target mass
$g$ will run, but exhibits significant tangling for very complicated images
like the world map, see Figure 8

Figure 8. Entire spherical mesh transformed by Optimal Information Transport 7
in 3D view. Shown on the right is detail on Africa, the Middle East, and Asia
with other side of sphere hidden from view in order to show the mesh more
clearly 8. One can observe some tangling here.

One potential solution to this issue exhibited by the Optimal Transport
problem will be explored in further work where the computation of the mapping
$T$ is slightly modified to gain convergence guarantees and greater stability.

###  5.2. Adaptive Mesh Methods for non-Spherical Blobs

Suppose that we have a mesh generated on an oblate sphere (ellipsoid of
revolution with eccentricity $e=0.25$), see 9

Figure 9. A $N=5048$-point cube mesh for the ellipsoid from the side 9 and
diagonally above looking down 9

For the Optimal Transport problem, the result from [9] asserts that there
exist $C^{\infty}$ source and target mass density functions $f_{0}$ and
$f_{1}$ such that the mapping $T$ is not guaranteed to be even continuous.
Further difficulty comes again from the geometry regarding the computation of
geodesics. Given a local coordinate chart $(x^{k})$ in a neighborhood of a
point $x$, we can solve the second order ordinary differential equation for
the geodesics:

(50)
$\ddot{x}^{k}(t)+\dot{x}^{i}(t)\dot{x}^{j}(t)\Gamma^{k}_{ij}\left(x(t)\right)=0$

where $\Gamma^{k}_{ij}$ are the Christoffel symbols. As the scheme stands,
however, we have exploited an explicit representation formula for the mixed
Hessian term. This required an explicit formula for the exponential map as a
surjective map from local tangent planes onto the sphere. This is no longer
possible in more general geometries. Thus, a direct discretization of
$D^{2}_{xy}c(x,y)$ would be necessary. Furthermore, the computation of the
mapping $T(x)$ can no longer depend on its exponential map representation.
That is, we must solve $\nabla u(x)=\nabla_{x}c(x,y)|_{y=T(x)}$. That is,
given $\nabla u(x)$ and $\nabla_{x}c(x,y)$, we find $T(x)$. We, similarly, do
not have a representation formula for geodesic normal coordinates. This means
that they must be approximated. This is possible with any method of finding
geodesic normal coordinates that is consistent.

The only difficulties that exist with the method described for Optimal
Information Transport are those concerning local geodesics. That means that
this method is more easily generalizable than that of Optimal Transport. These
observations will be further explored in future work.

##  6\. Conclusion

We have provided the first comparison of Optimal Transport and Optimal
Information Transport on the sphere for the diffeomorphic density matching
problem with applications to adaptive mesh methods by using convergent finite-
difference schemes. We also introduced the issues concerning their
generalization to more general closed, compact surfaces $M$. In this more
general case, there appear to be more merits to using Optimal Information
Transport. These merits include speed of implementation, generalizability in
computational terms but also in terms of a priori regularity results.
Nevertheless, it could be that a particular application requires the
computation of Optimal Transport with the squared geodesic cost on the sphere
and more general surfaces.

As an additional benefit, the numerical techniques presented here are more
computationally efficient than direct computation via discretization of the
probability measures via Dirac measures. Many of the numerical techniques
presented here apply more generally to elliptic PDE on compact closed
surfaces, potentially solving a wide variety of problems.

## References

  * [1] Martin Bauer, Sarang Joshi, and Klas Modin.  Diffeomorphic density matching by optimal information transport.  SIAM Journal on Imaging Sciences, 8(3):1718–1751, 2015.
  * [2] Martin Bauer, Sarang Joshi, and Klas Modin.  Diffeomorphic random sampling using optimal information transport.  Geometric Science of Information 2017, pages 135–142, 2017.
  * [3] C. J. Budd, M. J. P. Cullen, and E. J. Walsh.  Monge-Ampère based moving mesh methods for numerical weather prediction, with applications to the Eady problem.  Journal of Computational Physics, 236:247–270, 2013.
  * [4] C. J. Budd and JF Williams.  Moving mesh generation using the parabolic Monge-Ampère equation.  SIAM Journal on Scientific Computing, 31(5):3438–3465, 2009.
  * [5] Xavier Cabré.  Topics in regularity and qualitative properties of solutions on nonlinear elliptic equations.  Discrete and Continuous Dynamical Systems, 8(2):331–359, April 2002\.
  * [6] L. Chacón, G. L. Delzanno, and J. M. Finn.  Robust, multidimensional mesh-motion based on Monge-Kantorovich equidistribution.  Journal of Computational Physics, 230:87–103, 2011.
  * [7] Chong Chen and Öktem Ozan.  Indirect image registration with large diffeomorphic deformations.  SIAM Journal on Imaging Sciences, 11(1), 2017.
  * [8] M. G. Crandall, H. Ishii, and P.-L. Lions.  User’s guide to viscosity solutions of second order partial differential equations.  Bulletin of the American Mathematical Society, 27(1):1–67, July 1992.
  * [9] A. Figalli, L. Rifford, and C. Villani.  On the Ma-Trudinger-Wang curvature on surfaces.  Calculus of Variations, 39:307–332, 2010.
  * [10] B. D. Froese.  Meshfree finite difference approximations for functions of the eigenvalues of the Hessian.  Numer. Math., 138(1):75–99, 2018.
  * [11] David Gilbarg and Neil S. Trudinger.  Elliptic Partial Differential Equations of Second Order.  Springer, 2001.
  * [12] V. Gorbunova, J. Sporring, P. Lo, M. Loeve, H. A. Tiddens, M. Nielson, A. Dirksen, and M. de Bruijne.  Mass preserving image registration for lung CT.  Medical Image Analysis, 16:786–795, 2012.
  * [13] S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent.  Optimal mass transport for registration and warping.  Int. J. Comput. Vision, 60(3):225–240, 2004.
  * [14] B. F. Hamfeldt and A. G. R. Turnquist.  A convergence framework for optimal transport on the sphere.  arXiv preprint arXiv:2103.05739, 2021.
  * [15] Brittany Froese Hamfeldt and Axel G. R. Turnquist.  A convergent finite difference method for optimal transport on the sphere.  Journal of Computational Physics, 445, November 2021.
  * [16] Brittany Froese Hamfeldt and Axel G. R. Turnquist.  A convergent numerical method for the reflector antenna problem via optimal transport on the sphere.  arXiv preprint arXiv:2108.00088, July 2021.
  * [17] Hitoshi Ishii.  On the equivalence of two notions of weak solutions, viscosity solutions and distribution solutions.  Funkcialaj Ekvacioj, 38:101–120, 1995.
  * [18] G. Loeper.  On the regularity of solutions of optimal transportation problems.  Acta Mathematica, 202:241–283, 2009.
  * [19] G. Loeper.  Regularity of optimal maps on the sphere: The quadratic cost and the reflector antenna.  Archive for rational mechanics and analysis, 199(1):269–289, 2011\.
  * [20] G. Loeper and C. Villani.  Regularity of optimal transport in curved geometry: the nonfocal case.  Duke Math. J., 151(3):431–485, February 2010.
  * [21] X.-N. Ma, N. X. Trudinger, and X.-J. Wang.  Regularity of potential functions of the optimal transportation problem.  Archive for rational mechanics and analysis, 177(2):151–183, 2005\.
  * [22] R. J. McCann.  Polar factorization of maps on Riemannian manifolds.  Geom. Funct. Anal., 11:589–608, 2001.
  * [23] Klas Modin.  Generalised Hunter-Saxton equations, optimal information transport, and factorisation of diffeomorphisms.  The Journal of Geometric Analysis, 25:1306–1334, 2015.
  * [24] T A E Moselhy and Y M Marzouk.  Bayesian inference with optimal maps.  Journal of Computational Physics, 231:7815–7850, 2012.
  * [25] J. Moser.  One the volume elements on a manifold.  Trans. Amer. Math. Soc., 120:286–294, 1965.
  * [26] A. M. Oberman.  Convergent difference schemes for degenerate elliptic and parabolic equations: Hamilton–Jacobi equations and free boundary problems.  SIAM J. Numer. Anal., 44(2):879–895, 2006.
  * [27] Caleb Rottman, Martin Bauer, Klas Model, and Sarang C Joshi.  Weighted diffeomorphic density matching with applications to thoracic image registration.  5th MICCAI workshop on mathematical foundations of computational anatomy, 2015.
  * [28] C. Villani.  Topics in Optimal Transportation.  American Mathematical Society, 2003.
  * [29] C. Villani.  Stability of a $4$th-order curvature condition arising in optimal transport theory.  Journal of Functional Analysis, 255:2683–2708, 2008.
  * [30] Cédric Villani.  Optimal Transport: Old and New, volume 338 of A Series of Comprehensive Studies in Mathematics.  Springer, 2009.
  * [31] Hilary Weller, Philip Browne, Chris Budd, and Mike Cullen.  Mesh generation on the sphere using Optimal Transport and the numerical solution of a Monge-Ampère type equation.  Journal of Computational Physics, 308:102–123, March 2016.
  * [32] Medha V Wyawahare, Pradeep M Patil, and Hemant K Abhyankar.  Image registration techniques: An overview.  International Journal of Signal Processing, Image Processing and Patter Recognition, 2(3), September 2009.
  * [33] L. Younes.  Shapes and Diffeomorphisms, volume 171 of Applied Mathematical Sciences.  Springer-Verlag, Berlin, 2010.

## Appendix A: Modified Poisson

Recall the definition of the viscosity solution 4.

######  Theorem 4.

For any sequence $\epsilon_{n}\rightarrow 0$, $\epsilon_{n}>0$ and for $R>0$
large enough, the viscosity solution $w_{\epsilon_{n}}$ of the PDE

(51) $\max\left\\{-\Delta w_{\epsilon_{n}}(x)+f(x),\left\|\nabla
w_{\epsilon_{n}}(x)\right\|-R\right\\}+\epsilon_{n}w_{\epsilon_{n}}(x)=0,$

and for continuous $f(x)$ satisfying $\int_{M}f(x)dx=0$, converges uniformly
to the unique mean-zero Lipschitz solution $u$ of the PDE $-\Delta
u(x)+f(x)=0$.

###### Proof.

The equation (51) is proper, and therefore even though it is degenerate
elliptic, there exists a comparison principle, provided that $f(x)$ is
continuous, by [8]. By Perron’s method for continuous viscosity solutions,
there exists a unique solution $w_{\epsilon}$ to this equation for any
$\epsilon>0$.

The viscosity solution $u_{\epsilon}$ of the PDE $-\Delta
u_{\epsilon}+f(x)+\epsilon u_{\epsilon}=0$ is unique as well and is Lipschitz
for all $\epsilon>0$ provided that it remains uniformly bounded in $\epsilon$.
It is mean zero, that is $\int_{M}u_{\epsilon}=0$. The distributional solution
$\tilde{u}_{\epsilon}$ of
$-\Delta\tilde{u}_{\epsilon}+f(x)+\epsilon\tilde{u}_{\epsilon}=0$ is the same
as the viscosity solution $u_{\epsilon}=\tilde{u}_{\epsilon}$, see [17]. Thus,
in the sense of distributions, for any $v\in C^{2}$ and use the integration by
parts formula on the closed compact manifold $M$:

(52) $\int_{M}v\Delta u_{\epsilon}dx-\int_{M}u_{\epsilon}\Delta vdx=0$

Choosing $v\equiv 1$, we get:

(53) $\int_{M}\Delta u_{\epsilon}dx=0$

Hence, we integrate the equation where $u_{\epsilon}$ is known a priori to be
continuous:

(54) $\int_{M}\Delta u_{\epsilon}+f(x)+\epsilon u_{\epsilon}dx=0$

Thus, since $f$ is already chosen to integrate to zero, we get
$\int_{M}u_{\epsilon}dx=0$.

Fix $\epsilon>0$. We will use interior regularity estimates to show that the
Lipschitz constant of $u_{\epsilon}$ is uniformly bounded. The results in
Euclidean space can be used for the compact surface $M$ by using a
$C^{\infty}$-atlas to cover the manifold and applying the Euclidean interior
estimates in coordinate patches $\Omega_{i}$ that overlap by a fixed distance
$\delta>0$ if necessary. Expressed over small enough coordinate patches in
local coordinates, we use estimates that are valid for second-order uniformly
elliptic linear PDE with bounded first- and zeroth-order coefficients. We use
these estimates, because the Poisson equation expressed in a small enough
coordinate patch in local coordinates has first-order terms due to the
Christoffel symbols. However, since the coordinate patch is small enough,
these Christoffel symbols are bounded. Interior Schauder estimates of
uniformly elliptic linear PDE in non-divergence form with bounded coefficients
can be found in Corollary 6.3 of Gilbarg and Trudinger [11]:

(55)
$d\left\|Du_{\epsilon}\right\|_{C^{0}(\Omega^{\prime}_{i})}+d^{2}\left\|D^{2}u_{\epsilon}\right\|_{C^{0}(\Omega^{\prime}_{i})}\leq
C\left(\left\|u_{\epsilon}\right\|_{C^{0}(\Omega_{i})}+\left\|f\right\|_{C^{0}(\Omega_{i})}\right)$

where $d\leq\text{dist}(\Omega^{\prime},\partial\Omega)$. In order to avoid
the effect of the terms $d$ on the estimate, the result over the compact
manifold $M$ will be effected by choosing overlapping coordinate patches
$\Omega_{i}$ that have a uniform constant of overlap $\delta>0$. That is,
$d(x,y)>\delta$ for any $x\in M\setminus\Omega_{i},y\in M\setminus\Omega_{j}$
for $i\neq j$. Thus, we will be able to choose $d\geq\delta$. Doing this, we
derive an estimate over the whole compact manifold $M$:

(56)
$\delta\left\|Du_{\epsilon}\right\|_{C^{0}(M)}+\delta^{2}\left\|D^{2}u_{\epsilon}\right\|_{C^{0}(M)}\leq
C\left(\left\|u_{\epsilon}\right\|_{C^{0}(M)}+\left\|f\right\|_{C^{0}(M)}\right)$

The supremum on the right-hand side of (55) can be bounded by using the
Krylov-Safonov Harnack inequality in Euclidean space [5]. Denote
$v_{\epsilon}=u_{\epsilon}-\inf_{M}u_{\epsilon}$. Then, we have:

(57) $\sup_{B_{R}}v_{\epsilon}\leq
C\left(\inf_{B_{R}}v_{\epsilon}+\left\|f\right\|_{L^{2}(B_{2R})}\right)$

Covering the compact manifold by a finite covering of balls that of radius $R$
small enough so that the Christoffel symbols in these balls are bounded, we
get:

(58) $\sup_{M}v_{\epsilon}\leq C\left\|f\right\|_{L^{2}(M)}$

since we apply the fact that

(59) $\sup_{B_{1}\cup B_{2}}u_{\epsilon}\leq
C\left(\sup_{B_{1}}u_{\epsilon}+\left\|f\right\|_{L^{2}(M)}\right)\leq\tilde{C}\left(\inf_{B_{1}}u_{\epsilon}+\left\|f\right\|_{L^{2}(M)}\right)$

a finite number of times. Thus, since $u_{\epsilon}$ is mean-zero, we have
that $\left\|u_{\epsilon}\right\|_{C^{2,0}(M)}\leq C$ for a universal constant
$C$ that does not depend on $\epsilon$. Hence, we know that since the
$u_{\epsilon}$ is bounded in its derivative, in fact, the solutions coincide,
i.e. $w_{\epsilon}=u_{\epsilon}$, for $R$ chosen large enough.

By the stability of viscosity solutions [8], that
$\overline{U}=\limsup_{\epsilon}u_{\epsilon}$ is a viscosity subsolution of
$-\Delta u+f(x)$. Likewise, we can define $\underline{U}=\liminf u_{\epsilon}$
and it is a viscosity supersolution of $-\Delta u+f(x)$.

Take any sequence $\epsilon_{n}\rightarrow 0$ such that
$u_{\epsilon}\rightarrow U$ uniformly. Then, $U$ is the mean-zero viscosity
solution of $-\Delta u+f(x)$. For any sequence $\epsilon_{n}$, there exists a
subsequence $\epsilon_{n_{k}}$ such that $u_{\epsilon_{n_{k}}}\rightarrow U$
uniformly, by Arzela-Ascoli. Since every sequence contains a subsequence that
converges to $U$, we conclude that, in fact $u_{\epsilon_{n}}\rightarrow U$
for any sequence. ∎


Copyright for this paper by its authors. Use permitted under Creative Commons
License Attribution 4.0 International (CC BY 4.0).

FIRE’23: Forum for Information Retrieval Evaluation , December 15–18, 2023,
Goa University, Panjim, India

[orcid=0000-0002-0051-1535<EMAIL_ADDRESS>[orcid=0000-0001-7347-7086<EMAIL_ADDRESS>[orcid=0000-0002-0548-5751,
<EMAIL_ADDRESS>[orcid=0000-0003-3866-3013<EMAIL_ADDRESS>
# Extending Translate-Train for ColBERT-X to African Language CLIR

Eugene Yang  Dawn J. Lawrie  Paul McNamee  James Mayfield  Human Language
Technology Center of Excellence, Johns Hopkins University, Baltimore,
Maryland, USA

(2023)

###### Abstract

This paper describes the submission runs from the HLTCOE team at the CIRAL
CLIR tasks for African languages at FIRE 2023. Our submissions use machine
translation models to translate the documents and the training passages, and
ColBERT-X as the retrieval model. Additionally, we present a set of unofficial
runs that use an alternative training procedure with a similar training
setting.

###### keywords:

ColBERT-X Translate-Train PLAID CLIR JH POLO

##  1 Introduction

The HLTCOE team participated in all four CLIR tasks in the CIRAL shared task.
These CLIR tasks use English queries to search for Hausa, Somali, Swahili, and
Yoruba documents. Our systems primarily use PLAID [1], an implementation of
ColBERT [2] retrieval architecture that encodes each token as a vector. Prior
work has demonstrated successes in augmenting training queries and passages
with translation to match the CLIR target languages for training CLIR dense
retrieval models [3]. This technique, named Translate-Train, was evaluated in
widely spoken languages, such as Chinese, Persian, and Italian, with decent
machine translation models. In the recent 2022 TREC NeuCLIR track [4], a
ColBERT model called ColBERT-X, trained with Translate-Train, is the most
effective end-to-end neural dense retrieval model.111The ColBERT-X runs were
contributed by the organizers, and, thus, marked as manual runs. Although
unlikely, the performance might have been affected by knowledge that is only
accessible by the organizers.

This paper documents our attempt to adapt the Translate-Train training
technique in African languages, where machine translation models have
generally struggled. As summarized in Table 1, we submitted seven runs for
each CLIR task, including two runs using machine-translated documents with
BM25 with RM3 (Run 1) and English ColBERT (Run 2) [5]. Since Yoruba was not
included in the pretraining of XLM-RoBERTa, we introduced an additional masked
language model (MLM) fine-tuning using the Afriberta Corpus
[6]222https://huggingface.co/datasets/castorini/afriberta-corpus to enhance
the language model for the CIRAL tasks.

We use MS MARCO [7] training triples with English queries and machine-
translated African language passages (Hausa, Somali, Swahili, and Yoruba) to
perform Translate-Train for ColBERT-X based on the out-of-box XLM-RoBERTa-
Large model [8] (Run 4) and our MLM-fine-tuned version (Run 6). We also
compare with the ColBERT-X model trained with English MS MARCO, i.e., English-
Trained, (Runs 3 and 5) to understand the quality and usefulness of the
machine-translated MS MARCO in the African languages. Finally, we experiment
with a new technique, JH POLO [9], that uses large language models to generate
English training queries drawn from the retrieval collection to perform in-
domain retrieval fine-tuning (Run 7).

Table 1: Run Name Description.  | Run | Doc MT | Ret. Model | MLM FT | Ret. FT | JHPolo FT
---|---|---|---|---|---|---
(1) | dt.bm25-rm3 | ✓ | BM25+RM3 | — | — | —
(2) | dt.plaid | ✓ | ColBERT | ✗ | ET | ✗
(3) | plaid-xlmr.et | ✗ | ColBERT-X | ✗ | ET | ✗
(4) | plaid-xlmr.tt | ✗ | ColBERT-X | ✗ | TT | ✗
(5) | plaid-xlmr.mlmfine.et | ✗ | ColBERT-X | ✓ | ET | ✗
(6) | plaid-xlmr.mlmfine.tt | ✗ | ColBERT-X | ✓ | TT | ✗
(7) | plaid-xlmr.mlmfine.tt.jholo | ✗ | ColBERT-X | ✓ | TT | ✓

##  2 Machine Translation

We used automated machine translation (MT) in two principal ways for the
evaluation. First, we used document translation to create English-language
representations of the CIRAL document collections, as this directly enables
search using English queries. Second, we translated the MS MARCO passages from
English to the four African languages.

Transformer-based models were trained using Amazon’s Sockeye v2 toolkit [10]
with training data that was principally from the open source repository, OPUS
[11]. Preprocessing steps included: running the Moses tokenizer, removal of
duplicate lines, and learning of subword units using the subword-nmt toolkit.
Case was retained. Notable hyperparameters include: use of 6 layers in both
encoder and decoder; 512 dimensional embeddings; 8 attention heads; 2,048
hidden units per layer; 30,000 subword byte pair encoding (BPE) unit,
separately in source and target languages; batch size of 4,096; the Adam
optimizer with an initial learning rate of $2\times 10^{-4}$.

###  2.1 Document Translation

When the query language is known ahead of time it is possible to translate
documents into the query language, effectively reducing the CLIR problem to a
monolingual task. Of course the quality of automated machine translation can
vary considerably, and some queries can materially suffer if named-entities or
other essential query elements are mistranslated. When languages have fewer
resources, and when source and target languages differ in linguistic typology,
translation can be challenging.

To increase the likelihood of producing better quality document translations
we created synthetic training bitext so that our neural machine translation
models would have larger quantities of data to work with. In recently
published work McNamee and Duh [12] showed that back-translation can be
particularly efficacious in lower-resource settings, and helps with lexical
coverage in the resulting translation system. For this setup we first trained
English-to-Other models, and used these initial four models to back-translate
7 million sentences of web-crawled English news. Then for each language these
7 million synthetic translations were added to our human produced training
data (i.e., bitext from OPUS) to then train the forward models which were used
to create English language translations of the four African document
collections.

###  2.2 Translating MS MARCO

MS MARCO was created to support neural IR over English texts. To support the
Translate-Train approach for the cross language setting, we wanted to produce
translations of MS MARCO into Hausa, Somali, Swahili, and Yoruba. The original
English dataset consists of 8,841,823 passages containing 497 million words.
Table 2 shows the quantity of training bitext, translation quality scores on a
commonly used benchmark, and the number of words in the translated MS MARCO
dataset, by language.

Table 2: MS MARCO Translations. Shown are: (a) the size of training bitext, in sentences; (b) MT quality using BLEU scores (lower-cased sacrebleu) on the FLORES-101 test set; and, (c) the size in words (millions) of the resulting translation. Language | Bitext size | FLORES-101 | Translation size
---|---|---|---
Hausa | 2.2M | 26.1 | 576M
Somali | 786k | 13.6 | 559M
Swahili | 9.9M | 37.7 | 502M
Yoruba | 1.4M | 5.5 | 672M

##  3 Training Pipeline

Our full training pipeline for ColBERT-X starts from the pretrained XLM-
RoBERTa Large model, followed by masked language model fine-tuning (MLM),
retrieval fine-tuning with translate-train, and finally, in-domain fine-tuning
with JH POLO. This section describes each fine-tuning step.

###  3.1 Masked Language Model Fine-tuning

Since XLM-RoBERTa [8] pretraining does not include Yoruba, we designed a fine-
tuning step to accommodate this absence. However, presenting only Yoruba text
to the model during fine-tuning risks catastrophic forgetting of other
language knowledge. Specifically, we would like the language model to retain
language knowledge related to the four African languages and to the query
language – English. Therefore, we present documents in Hausa, Somali, Swahili,
Yoruba, and English round-robin to perform masked language model fine-tuning.
We used Common Crawl documents in Afriberta Corpus [6] for the four African
languages and collected additional English Common Crawl documents to match the
genre.

We fine-tune the model for 200,000 update steps using a learning rate of
$1\times 10^{-5}$ and a batch size of 48 text sequences of a maximum length of
512 tokens each. We used four A100 NVidia GPUs to train the model. Fine-tuning
took around 34 hours to complete.

###  3.2 Retrieval Fine-tuning with Translate-Train

To transform a multilingual language model into a CLIR ColBERT-X model, we
fine-tuned the language model using MS MARCO small training triples with the
original English queries and translated passages (Translate-Train) [3]. We
evaluate this Translate-Train with both the pretrained XLM-RoBERTa model and
our MLM-fine-tuned language model. The model is trained with a contrastive
loss using Cross-Entropy between the positive and negative passages of each
training query. For comparison, we also fine-tuned the language model with
English MS MARCO without translation (English-Train).

We fine-tune the language model with the retrieval objective for 200,000
update steps with a learning rate of $5\times 10^{-6}$ and a batch size of 64
triples (query, positive, and negative passage triplets). Following the
ColBERT-X [3] training setup, we pad the queries to 32 tokens with [MASK]
tokens. Each ColBERT-X model is trained with eight V100 NVidia GPUs for around
50 hours. For the official submissions, we used the PLAID [1] implementation
of ColBERT training. However, after the submission, we discovered that the
ColBERT-X implementation 333https://github.com/hltcoe/ColBERT-X, which is
based on the ColBERT v1 [2] codebase, provides a more stable and effective
training process. Thus, we also report a set of unofficial runs using this
implementation.

For retrieval, we use the PLAID [1] retrieval implementation, which uses
K-Means clustering and compression to approximate and accelerate retrieval. We
compress each document token residual vector dimension down to one bit,
resulting in a 128-bit residual representation for each document token.

###  3.3 JH POLO In-Domain Retrieval Fine-tuning

Training data for the CIRAL languages is quite limited. One option for new
training data is Translate-Train: translating the documents of an existing
retrieval training collection, such as MS MARCO, to the target languages.
However, machine translation for the CIRAL languages is not particularly good
at the time of the evaluation. Furthermore, there is no guarantee that the
documents of an existing evaluation collection will be a good match for those
of the target collection. Creating new training examples using the target
collection itself for the documents would eliminate these problems; documents
would be naturally occurring, and would therefore not exhibit
“translationese.” And there would never be a mismatch between the genre or
style of the documents in the training collection and those in the target
collection.

JH POLO [9] is a methodology for creating such training data. It relies on the
existence of a large generative language model that includes coverage of the
target language. The process begins by selecting two documents from the target
collection that have some topic overlap. One of the documents will end up as a
relevant document in a training example, and the other will become a non-
relevant document in the same example. Selecting document pairs that are
closer in meaning will lead to harder negative examples in the training
examples produced.

Once the documents have been selected, the generative language model is
prompted to create a query for which the first document is relevant and the
second document is not. This query, and the two documents, are bundled to form
a single new training example. This process can be repeated to generate as
many training examples as desired.

You must write questions for a news quiz to appear in the newspaper. A news
quiz asks about events in the news, NOT about news articles. Here are two
articles that appeared in this week’s news: <<first>> <<second>> For each
article give five factual news quiz English questions, one per line with no
extraneous words, that are answered by the events described in that document
and are not answered by the events described in the other document. The quiz
questions must never refer to individual news articles, or assume the quiz-
taker has seen those articles. Precede the first five with DOCA: and the
second with DOCB: Figure 1: GPT-4 prompt used to create JH POLO training
examples.

We used the JH POLO methodology to create training data for the four CIRAL
languages. We used GPT-4 for the generative language model. While GPT-4 would
occasionally complain that it was unable to handle documents in one of the
CIRAL languages, in almost all cases it would willingly process the documents
without being told what language they were written in. In addition to allowing
naturally-occurring documents in the training set, this approach hits the
sweet spot of most generative LLMs: producing short, English texts. Our prompt
is shown in Figure 1. The prompt accomplishes several things:

  * •

It identifies the task as question answering.

  * •

It attempts to focus the questions on the content of the news articles, not on
the articles themselves. It also tries to prevent the generated queries from
assuming that the searcher knows the content of any article a priori. GPT-4
had a difficult time conforming to these requirements.

  * •

It includes the text of the two documents or passages.

  * •

It asks for ten training examples for each document pair, with each document
serving as the relevant document for five of thOSe queries.

  * •

It specifies how the output should be formatted.

⬇ {  "id": "chatcmpl-7ruApne33HFrrRiJB25563Ctcr3Ud",  "object":
"chat.completion",  "created": 1693081671,  "model": "gpt-4-0613",  "choices":
[  {  "index": 0,  "message": {  "role": "assistant",  "content": "DOCA:\n1.
Who reflected on President Buhari’s leadership in the recent public
discourse?\n2. What is the criticism regarding the opposition’s treatment
raised in the recent debate?\n3. What were the financial implications
mentioned in recent political critiques?\n4. What are the societal issues
addressed by the speaker?\n5. Who made the prayer for wisdom and understanding
at the end of the recent speech?\n\nDOCB:\n1. What were the tragic aspects of
Yunusa and Ese’s love story mentioned in recent reports?\n2. Who recently
criticized some religious leaders for committing wrongdoings?\n3. What legal
judgement was recently confirmed as punishment for an offender?\n4. What issue
of child exploitation came to light recently?\n5. What phrase has been adopted
by vocal sympathizers to describe the prevailing situation?"  },
"finish_reason": "stop"  }  ],  "usage": {  "prompt_tokens": 742,
"completion_tokens": 161,  "total_tokens": 903  }  }  Figure 2: GPT-4 output
used to create JH POLO training examples.

To select document pairs, we first used each document with more than a fixed
number of characters (the “query document”) as a query over the other
documents using a BM25 sparse retrieval model. We considered each of the top
twenty documents in the resulting ranked list (the “candidate document”) not
including the query document. We eliminated from consideration any document
that met any of the following criteria:

  * •

the ratio of the score of the candidate document to that of the query document
was greater than 0.65

  * •

the longest common substring between the query document and the candidate
document was more than 60% of the entire candidate document

  * •

fewer than twenty characters from the candidate document were not part of the
longest common substring

  * •

the candidate document had fewer than 150 characters

We selected for inclusion in the training collection the pair that was not
rejected by the above criteria, and that maximized the size of the training
collection, given that no document was allowed to be part of more than one
pair.

Once the document pairs were selected, we embedded the text of each document
in the GPT-4 prompt and ran the prompt. In most cases, GPT-4 successfully
produced output with ten output queries per prompt. Figure 2 shows the GPT-4
output for a completed prompt.

We applied two forms of automated quality control to the JH POLO outputs.
First, because GPT-4 had a difficult time omitting mention of the documents in
its output queries and not assuming the user knew anything about those
documents, we eliminated any query that contained any of the words articles,
reports, speaker, these. Second, to try to eliminate examples where the
relevant and non-relevant documents were too close together, we used an
mMiniLM cross-encoder (cross-encoder/mmarco-mMiniLMv2-L12-H384-v1) to compare
the query to each of the documents; we eliminated any example where the cross-
encoder score (between 0 and 1) for the positive document was not at least
0.15 above the score of the non-relevant document. The result was a collection
of 48,459 training examples over 14,323 document pairs in the four CIRAL
languages combined.

##  4 Results

Table 3 summarizes the run results of our official submissions. For all four
languages, an English ColBERT model indexing the machine-translated documents
(Run 2) provides the most effective retrieval results in both nDCG@20 and
R@100. For Somali, this model ranks the best results among all submissions.

The ColBERT-X models submitted are not as effective as the monolingual English
one. Similar to prior works, models trained with Translate-Train are more
effective than those trained with English-Train. The difference in Yoruba
without MLM fine-tuning is substantially larger (Runs 3 and 4), since
Translate-Train is the only step that introduces Yoruba text to the model
during training.

Table 3: Official runs from the HLTCOE team. Runs are evaluated by the track
organizers.

|  |  |  | nDCG@20 | R@100
---|---|---|---|---|---
| MLM | Ret. FT | JH POLO | Hausa | Somali | Swahili | Yoruba | Hausa | Somali | Swahili | Yoruba
| Max | 0.5700 | 0.5118 | 0.5232 | 0.5819 | 0.5902 | 0.6436 | 0.5956 | 0.8057
| Median | 0.2530 | 0.2445 | 0.2447 | 0.3138 | 0.3576 | 0.3083 | 0.3340 | 0.5037
| Mean | 0.2690 | 0.2403 | 0.2644 | 0.3115 | 0.3598 | 0.3265 | 0.3249 | 0.5091
(1) | DT >> BM25+RM3 | 0.2015 | 0.2550 | 0.2178 | 0.3555 | 0.3359 | 0.4210 | 0.3340 | 0.6273
(2) | DT >> English ColBERT | 0.4743 | 0.5118 | 0.4932 | 0.4793 | 0.5733 | 0.6436 | 0.5956 | 0.7240
ColBERT-X Submission Models |
(3) | ✗ | ET | ✗ | 0.3481 | 0.2915 | 0.3182 | 0.2627 | 0.5237 | 0.4332 | 0.4616 | 0.5176
(4) | ✗ | TT | ✗ | 0.3557 | 0.2878 | 0.3347 | 0.3522 | 0.5083 | 0.4373 | 0.4510 | 0.5784
(5) | ✓ | ET | ✗ | 0.3488 | 0.2760 | 0.3301 | 0.3804 | 0.5326 | 0.4260 | 0.4487 | 0.6950
(6) | ✓ | TT | ✗ | 0.4335 | 0.3366 | 0.4230 | 0.4189 | 0.5256 | 0.4534 | 0.4645 | 0.6394
(7) | ✓ | TT | ✓ | 0.3601 | 0.3117 | 0.4081 | 0.4297 | 0.4829 | 0.4277 | 0.4477 | 0.6748

We observe that MLM fine-tuning is generally helpful even when followed by
Translate-Train, which also conveys African language knowledge to the model.
The MLM-fine-tuned model followed by English-Train demonstrates similar
effectiveness when directly performing Translate-Train on the out-of-box XLM-
RoBERTa Large model. This observation aligns with prior work in continued
language model fine-tuning for CLIR tasks [13], where the authors found that
an effective language model fine-tuning step can replace Translate-Train.
Furthermore, performing Translate-Training on top of the MLM-fine-tuned
language model also leads to better effectiveness, suggesting THAT all
training steps contribute to the final retrieval effectiveness.

However, the in-domain JH POLO fine-tuning does not seem to be helpful. For
Hausa, Somali, and Swahili, fine-tuning on an additional 1,000 JH POLO
training examples degraded performance. However, both nDCG@20 and R@100
improved after JH POLO fine-tuning for Yoruba. We hypothesize that this is due
to the additional Yoruba text that is presented to the model. Since the model
has only seen a small amount of Yoruba text, any additional training showing
more Yoruba language would be helpful.

###  4.1 Unofficial Runs

Table 4: Unofficial Runs. Runs are evaluated with the qrels provided by the
organizers and evaluated by the HLTCOE team.

|  |  | nDCG@20 | Judged@20
---|---|---|---|---
MLM | Ret. FT | JHPolo | Hausa | Somali | Swahili | Yoruba | Hausa | Somali | Swahili | Yoruba
✗ | ET | ✗ | 0.3516 | 0.3080 | 0.3064 | 0.2534 | 0.5225 | 0.4763 | 0.5271 | 0.4995
✗ | TT | ✗ | 0.3722 | 0.2994 | 0.3766 | 0.3625 | 0.5925 | 0.5354 | 0.6165 | 0.6490
✓ | ET | ✗ | 0.3751 | 0.3018 | 0.3179 | 0.4097 | 0.5563 | 0.4848 | 0.5382 | 0.6475
✓ | TT | ✗ | 0.3450 | 0.2513 | 0.3093 | 0.3863 | 0.5562 | 0.4515 | 0.5124 | 0.6535
✓ | ET | ✓ | 0.3451 | 0.3069 | 0.3083 | 0.4105 | 0.5369 | 0.4914 | 0.5147 | 0.6665
✓ | TT | ✓ | 0.2957 | 0.2276 | 0.2878 | 0.4168 | 0.4406 | 0.3884 | 0.4682 | 0.6400

Based on other experiments, we discovered that the PLAID training
implementation (essentially version 3 of the ColBERT implementation) leads to
degraded performance in the resulting IR model. We retrained the models using
the original ColBERT-X implementation and present the results in Table 4.
Since these runs are produced after the submission deadline, the runs are not
part of the pooling assessments. Therefore, only around 50% to 60% of the top
20 retrieved documents are judged. While treating the unjudged documents as
non-relevant is a common assumption in IR evaluation, this also suggests that
results presented in Tables 3 (and other official submissions) and 4 are not
perfectly comparable.

Based on results in Table 4, models trained with the ColBERT-X implementation
seem to be generally more effective. While the trend of the contribution
provided by each training step is less clear, Translate-Train without MLM
still provides more effective models than English-Train, except for Somali.

However, based on this set of results, the benefit of the additional MLM fine-
tuning step is smaller. In fact, the knowledge in the Afriberta Corpus and in
the machine-translated MS MARCO seem to be contradictory. While performing
only Translate-Train or MLM fine-tuning still leads to similar effectiveness,
doing both does not give us additional advantage.

## References

  * Santhanam et al. [2022a] K. Santhanam, O. Khattab, C. Potts, M. Zaharia, PLAID: an efficient engine for late interaction retrieval, in: Proceedings of the 31st ACM International Conference on Information & Knowledge Management, 1747–1756, 2022a.
  * Khattab and Zaharia [2020] O. Khattab, M. Zaharia, Colbert: Efficient and effective passage search via contextualized late interaction over bert, in: Proceedings of the 43rd International ACM SIGIR conference on research and development in Information Retrieval, 39–48, 2020.
  * Nair et al. [2022] S. Nair, E. Yang, D. Lawrie, K. Duh, P. McNamee, K. Murray, J. Mayfield, D. W. Oard, Transfer Learning Approaches for Building Cross-Language Dense Retrieval Models, in: Advances in Information Retrieval: 44th European Conference on IR Research, ECIR 2022, Stavanger, Norway, April 10–14, 2022, Proceedings, Part I, Springer-Verlag, Berlin, Heidelberg, ISBN 978-3-030-99735-9, 382–396, URL https://doi.org/10.1007/978-3-030-99736-6_26, 2022\.
  * Lawrie et al. [2023] D. Lawrie, S. MacAvaney, J. Mayfield, P. McNamee, D. W. Oard, L. Soldaini, E. Yang, Overview of the TREC 2022 NeuCLIR Track, 2023.
  * Santhanam et al. [2022b] K. Santhanam, O. Khattab, J. Saad-Falcon, C. Potts, M. Zaharia, ColBERTv2: Effective and Efficient Retrieval via Lightweight Late Interaction, in: Proceedings of the 2022 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Association for Computational Linguistics, Seattle, United States, 3715–3734, URL https://aclanthology.org/2022.naacl-main.272, 2022b.
  * Ogueji et al. [2021] K. Ogueji, Y. Zhu, J. Lin, Small Data? No Problem! Exploring the Viability of Pretrained Multilingual Language Models for Low-resourced Languages, in: Proceedings of the 1st Workshop on Multilingual Representation Learning, Association for Computational Linguistics, Punta Cana, Dominican Republic, 116–126, URL https://aclanthology.org/2021.mrl-1.11, 2021.
  * Nguyen et al. [2016] T. Nguyen, M. Rosenberg, X. Song, J. Gao, S. Tiwary, R. Majumder, L. Deng, MS MARCO: A Human Generated MAchine Reading COmprehension Dataset, CoRR abs/1611.09268, URL http://arxiv.org/abs/1611.09268.
  * Conneau et al. [2020] A. Conneau, K. Khandelwal, N. Goyal, V. Chaudhary, G. Wenzek, F. Guzmán, E. Grave, M. Ott, L. Zettlemoyer, V. Stoyanov, Unsupervised Cross-lingual Representation Learning at Scale, in: Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, Association for Computational Linguistics, Online, 8440–8451, URL https://aclanthology.org/2020.acl-main.747, 2020\.
  * Mayfield et al. [2023] J. Mayfield, E. Yang, D. Lawrie, S. Barham, O. Weller, M. Mason, S. Nair, S. Miller, Synthetic Cross-language Information Retrieval Training Data, arXiv preprint arXiv:2305.00331 .
  * Vaswani et al. [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, I. Polosukhin, Attention is All You Need, in: Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, Curran Associates Inc., Red Hook, NY, USA, ISBN 9781510860964, 6000–6010, 2017\.
  * Tiedemann [2012] J. Tiedemann, Parallel Data, Tools and Interfaces in OPUS, in: Proceedings of the Eighth International Conference on Language Resources and Evaluation (LREC’12), European Language Resources Association (ELRA), Istanbul, Turkey, 2214–2218, URL http://www.lrec-conf.org/proceedings/lrec2012/pdf/463_Paper.pdf, 2012\.
  * McNamee and Duh [2023] P. McNamee, K. Duh, An Extensive Exploration of Back-Translation in 60 Languages, in: Findings of the Association for Computational Linguistics: ACL 2023, Association for Computational Linguistics, Toronto, Canada, 8166–8183, URL https://aclanthology.org/2023.findings-acl.518, 2023\.
  * Yang et al. [2022] E. Yang, S. Nair, R. Chandradevan, R. Iglesias-Flores, D. W. Oard, C3: Continued Pretraining with Contrastive Weak Supervision for Cross Language Ad-Hoc Retrieval, in: Proceedings of the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval, SIGIR ’22, Association for Computing Machinery, New York, NY, USA, ISBN 9781450387323, 2507–2512, URL https://doi.org/10.1145/3477495.3531886, 2022\.